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ABSTRACT 


A pool of liquid with a horizontal free surface is bounded on one side by a 
vertical solid wall, which is maintained at a cold temperature relative to the core flow 
region. Strong temperature gradients along the surface give rise to surface tension 
variations (thermocapillary stress), which drives flow. Thin viscous boundary layers 
form along the surface and wall. A boundary-layer model is designed which captures 
the dynamics of the cold corner, applicable for any Marangoni number M and Prandtl 
number P in the convective inertial regime. 

Analytical expressions for the velocity and boundary-layer thicknesses are de¬ 
veloped, which allow accurate prediction of the flow field. The core flow region (out¬ 
side the viscous boundary layers) is treated as irrotational flow and Laplace’s equation 
is solved using both a Green’s function approach and a complex variables approach 
in the quarter-plane. The flow along the wall is treated as a plane wall jet. 

The two-dimensional unsteady heat equation is solved using an alternating 
direction implicit method. Results show that the flow into the corner is strong enough 
to contain the thermal field, compressing the isotherms along the wall after steady- 
state is reached. Additionally, a uniform stream function prediction is developed, by 
matching the inner and outer flows, giving a relatively accurate depiction of the flow. 
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DISCLAIMER 


The computer programs in Appendix E are supplied on an “as is” basis, with 
no warrantees of any kind. The author bears no responsibility for any consequences 
of using these programs. 
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I. 


INTRODUCTION 


One must learn by doing the thing; for though you think you know it, you 

HAVE NO CERTAINTY UNTIL YOU TRY. — Sophocles 


The purpose of this research is to exa mi ne the velocity and temperature fields 
in a two-dimensional fluid flow problem embodying a thermocapillary feedback mech¬ 
anism. The structure of the velocity and thermal fields will change, based upon 
whether the heat transfer is conductive or convective in nature, as well as whether 
inertia plays a dominant role or not. One of the governing parameters, the Marangoni 
number M, measures the importance of thermal convection relative to thermal dif¬ 
fusion. A second governing parameter, the Prandtl number P, gives the ratio of 
viscous to thermal diffusion for the material. In this research, we study the effects 
of a high Marangoni number, indicative of convective heat transfer, combined with 
a low Prandtl number, indicative of inertial flow. The most important regime for 
materials processing is that regime where thermal convection is dominant and the 
flow is inertial, with viscous boundary layers for min g. This dissertation out lin es a 
model for this important regime, with the goal of predicting the velocity and thermal 
fields which will occur. 

Chapter II gives a background on the importance of this research and possible 
physical applications, followed by a survey of related work in thermocapillaxy flows by 
previous researchers. In our literature search, we did not find any other researchers 
who had studied the cold comer region with large Marangoni numbers and small 
Prandtl numbers. A representative velocity vector field (using Canright’s numerical 
data [Ref. 1]), which serves as a guiding principle for our assumptions, is shown at 
the chapter’s end. 

Chapter III formulates the problem mathematically and gives all assumptions 
and a thorough scaling analysis. It is important to obtain a summary of theoretical 
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scales of the thermal length, viscous thickness and length, and surface velocity. Al¬ 
though the scalings derived are different than those previously published, they agree 
better with the published numerical data. 

Chapter IV contains the conservation equations and how they affect each of 
four regions within the cold corner. Using the appropriate scaling factors, once the 
velocity is determined in one region, it is matched to that in all neighboring regions. 

Chapter V discusses the approximation methods employed to describe the 
velocity profiles in the viscous boundary layers and core flow region. Two methods 
from the literature were modified to predict velocity behavior along the free surface 
and the rigid wall. The modification along the free surface introduces a new technique 
to predict the velocity from a free surface with a thermal gradient boundary condition 
and no non-slip condition through a viscous boundary layer into the core flow region. 
The flow from the corner down the rigid wall is treated as a plane wall jet. An 
expression for the temperature solution in the domain is then determined. 

Chapter VI shows the numerical results obtained from various techniques. 
The alternating direction implicit (ADI) method is the basis for solving the two- 
dimensional unsteady heat equation. Various tests against problems with known 
solutions are performed to validate the code, followed by the solution of the present 
problem. The scheme is also modified to update the velocity profile during the time- 
step iterations to model the physical feedback. A modification to a wall jet profile is 
incorporated into the scheme. The results are presented as a uniform velocity field, 
combining the boundary layers and core flow, for comparison with previous numerical 
results. 

Chapter VII gives the conclusions and a few discussions associated with them. 
Chapter VIII fists selected areas of possible future research. 

Appendices A, C, and D show work which attempted to improve the model in 
the surface and wall boundary layers. We ultimately did not use these methods (for 
reasons given in each appendix), but we felt we should include them. Appendix B 
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gives a derivation of an important and often-used constant. The numerical codes are 
listed in Appendix E. 
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II. BACKGROUND 


We do not set out on this path of enquiry just as adventurers. We follow 

THE FOOTSTEPS OF SUCH GIANTS AS EULER AND LAGRANGE, WHO ALSO WONDERED 
WHERE IT LED, AND PRANDTL, WHO WALKED IT AND FOUND THE ANSWER. — David 
Pnueli and Chaim Gut finger 


A. MOTIVATION 

Materials processing often involves the melting and solidification of the mater¬ 
ial. Several practical processes, such as welding, containerless processing of materials, 
float-zone purification, and Czochralski crystal growth, involve a pool of molten metal 
with a free surface [Ref. 1]. The purest large-volume silicon single crystals are being 
commercially grown by floating-zone melting [Ref. 2, 3]. Inevitably, there exists a 
fluid flow within this molten region. 

When heavier fluid is on top of fighter fluid, the weight of the heavier, denser 
fluid makes it tend to sink, while the fighter, less dense fluid tends to rise. This sets up 
motion within the fluid known as buoyancy-driven convection or natural convection , 
also known as Rayleigh convection. One of the most commercially important materials 
processes that are subject to buoyancy-driven convection is that of crystal growth. 
A major advantage of performing experiments in a space laboratory is being able to 
minimize gravitational effects and therefore reduce buoyancy-driven convection [Ref. 
4]. There are important reasons to avoid having buoyancy-driven convection during 
crystal growth. Convection can alter trace element distribution in the resulting crystal 
and may degrade its quality. Many other physical processes can greatly affect crystal 
growth, such as radiation heat transfer and surface tension driven flow. Typically, all 
of these processes happen at the same time on Earth. Buoyancy-driven convection is 
the predominant fluid flow on Earth, since the gravitational forces are stronger than 
surface tension forces [Ref. 5]. In microgravity, however, buoyancy-driven convection 
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is reduced by a factor of 10 -6 , allowing thermocapillary flow to predominate [Ref. 
6]. The near elimination of buoyancy-driven convection has promise for materials 
processing in space, a low-gravity environment. For such applications, it is essential 
that thermocapillary flow be studied in detail, without the influence of buoyancy- 
driven convection. 

Convection in the molten metal is typically vigorous and can greatly affect 
the results of the materials process, including the size and shape of the melt pool, 
the heat transfer, the mixing of solutes, and the final microstructure of the product 
[Ref. 7]. The buoyancy force can also play a role in establishing the fluid flow. 
When the dimension of the molten region is small, as in the case of welding and 
laser materials processing, the buoyancy force becomes unimportant. Beyond the 
molten region, the absorbed heat is conducted into the substrate and transported 
away by the motion of the substrate. Surface tension gradient driven flow has been 
identified to be responsible for ripple formation in the weld [Ref. 7]. When buoyancy, 
electromagnetic, and surface tension forces were considered in a numerical solution 
for the convective heat transfer within the molten pool during arc welding, it was 
found that the surface tension gradient is the dominant factor in many cases [Ref. 
8]. Other authors have verified the dominance of thermocapillary flows in materials 
processing [Ref. 9, 10, 11]. 

In addition to the heat transport due to the motion of the substrate, there 
is also the convection due to the fluid flow within the molten pool. Depending on 
the type of heat source, the fluid flow may be caused by the surface tension gradient 
[Ref. 7]. When a free surface is heated by a concentrated heat source, the resulting 
temperature distribution causes a nonuniform surface tension distribution [Ref. 12], 
Typically, surface tension is a decreasing function of temperature. Thus, the fluid 
layer on the surface experiences a shear force pulhng the fluid from the warmer central 
area of the molten pool to the cooler outer region. Since the bulk fluids are viscous, 
they are dragged along; bulk fluid motion then results from the temperature gradients 
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along the interface. This phenomenon is known as thermocapillary or Marangoni 
convection. The Marangoni effect is associated with movement in a fluid interface, 
caused by local variations in interfacial tension that are caused in turn by differences 
in composition or temperature. A drying coat of paint is a simple example where local 
variations of surface tension set into intricate motion an entire liquid film [Ref. 13]. In 
the case of the float-zone crystal growth process, a melt zone is moved slowly through 
a cylindrical, vertically-oriented crystal. The liquid zone between the polycrystalline 
incoming material and the single crystal that is grown is stabilized by surface tension 
[Ref. 3]. The temperature distribution varies along the interface of the liquid and gas, 
which leads to surface-tension gradients, creating a considerable thermal Marangoni 
flow driving force, or thermocapillary convection. 

B. PREVIOUS RESEARCH 

There have been many theoretical studies of thermocapillary flows, both nu¬ 
merical and analytical. Ostrach [Ref. 5] reviewed of some of the analytical studies rel¬ 
evant to low-gravity conditions. Some fluid flows, such as buoyancy-driven and g-jitter 
convection, are acceleration dependent, while others, such as surface-gradient, ther¬ 
moacoustic, and phase-change convection are independent of the acceleration field. 
There are forces in a space laboratory which act on experiments in a gravity-like man¬ 
ner. G-jitter (or residual acceleration [Ref. 4]) is the term used when combining these 
gravity-like forces together. Thermoacoustic convection refers to the pressure-driven 
convection that results when a confined compressible fluid is heated rapidly; the nam e 
describes the sonic character of induced pressure waves. When a fluid solidifies or 
vaporizes, a change in density associated with the phase transition may also cause a 
change in volume. A shrinkage occurring during solidification may result in a volume 
reduction, further resulting in the flow of a liquid toward the solidifying surface. This 
is known as phase-change convection. Davis [Ref. 9] reviewed thermocapillary in¬ 
stabilities in planar layers. He gave details about studies concerning Marangoni and 
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hydrothermal instabilities involving the one-layer system geometry. Kuhlmann [Ref. 
3] gave a thorough reference list of many processes subject to Marangoni effects, such 
as the motion of droplets and bubbles, the stabilization of thin films, evaporation 
and boiling, welding, electrokinetic effects, and crystal growth techniques, to name 
just a few. Hondros [Ref. 14] offered several examples of phenomena attributed to 
Marangoni convection relevant to modern technology, such as in hot salt corrosion in 
turbine blades, the drying of silicon wafers in the electronics industry, and micro-pools 
produced by plasma disruptions in a prospective thermonuclear fusion reactor. 

Cowley and Davis [Ref. 15} analyzed the two-dimensional thermocapillary flow 
near a hot wall for vigorous, viscous flow (large Marangoni and Prandtl numbers). 
The fluid flows up a wall, which is kept at a constant high temperature within a finite 
distance from the free surface, and then turns and flows along the free surface. This 
is called the hot corner problem. The finite distance determines the length scale. The 
parameter regimes studied either have no viscous boundary layers or have viscous 
boundary layers that are much thicker than the thermal layers. The solution they 
find is the thermocapillary analog of buoyancy-driven convection in a quarter plane 
described by Roberts [Ref. 16]. 

Sen and Davis [Ref. 11] studied steady flow in two-dimensional slots, which are 
differentially heated, inducing the surface-tension gradient along the interface. They 
studied the case with the aspect ratio A (ratio of the depth to the width) approaching 
zero. The ends of the slots are maintained at a fixed temperature difference. Flows 
on the interfaces are directed from the hot towards the cold end and return along 
a region removed from the interfaces. The pressure is higher at the cold end and 
the interface thus bulges near the cold end and is constricted near the hot end. The 
leading-order outer solutions having parallel flow and flat interfaces continue to be the 
leading order outer approximations even when conditions on the Reynolds number 
and Marangoni number are relaxed from 0(A). 

Zebib, Homsy, and Meiburg [Ref. 17] conducted numerical studies of flow in 
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a two-dimensional square cavity with one hot wall and one cold wall. They showed 
that for moderate to small Prandtl numbers, the cold corner has the stronger effect 
on both the flow and the heat transfer of the cavity. Their results at finite Marangoni 
number in a finite pool suggest that experimentally observed instabilities may be 
associated with the rapid turning flows and high vorticity in the cold wall region. 
Another important finding in this paper was the determination of the scalings for the 
thermocapillary boundary-layer thickness on the free surface and on the rigid wall. 
Their results confirmed those of Ostrach [Ref. 5, 24]. 

Ohring and Lugt [Ref. 10] studied the onset of Marangoni convection in a float 
zone of liquid silicon from a state at rest in the absence of gravity, with high Marangoni 
numbers. They found the existence of a critical Marangoni number at which the 
steady flow field becomes unstable and changes to an axisymmetric oscillatory field. 
Also, they concluded that a deformable free surface, even a very small one on the 
order of 10~ 4 cm, is necessary for the onset of instability. The velocity at the free 
surface is extremely high. The temperature field oscillates and causes uneven heat 
transfer at the walls, which is not desirable for silicon crystal growth. 

Chan, Mazumder, and Chen [Ref. 18] developed a three-dimensional model of 
the fluid flow and heat transfer of a laser melted pool. Using a perturbation solution, 
they found that the presence of the thermocapillary convection causes the physics of 
the problem to change from conduction to convection dominated. The pool geometry 
then changes dramatically, resulting in up to a 150 percent increase in the aspect ratio 
as compared to the pure conduction case. In another paper [Ref. 12], they presented 
solutions of thermocapillary convection in the central region of a nonuniformly heated 
surface for the asymptotic limits of very high and very small Prandtl numbers. The 
model is valid when the viscous and thermal boundary layers are small compared to 
the depth and width of the melt pool. In addition, the intensity of the thermocapil¬ 
lary convection and the boundary-layer thicknesses in the stagnation region depend 
primarily on the curvature of the heat flux distribution. In a comprehensive report, 
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Mazumder, Chen, Chan, and Zehr [Ref. 7] developed a three-dimensional perturba¬ 
tion model for surface tension driven flow with a flat surface which predicted velocity 
field, temperature field, effect of trace elements and a self-consistent prediction of the 
pool shape and cooling rate. The three-dimensional problem is represented by two 
sets of two-dimensional governing equations. It is found that a particle recirculates 
many, many times in the molten pool before it resolidifies, showing that the solute 
can be well mixed. Regarding free surface deformation, thermocapillarity drives the 
surface fluid radially outward at extremely high velocities. These high velocities dis¬ 
place more mass from the surface region than can be replaced by the recirculating 
flow, which thus causes a depression. The displaced mass builds up at a solid/liquid 
interface, which causes the surface to bulge upward where the liquid turns downward 
into the molten pool. Experiments were performed to examine the validity of the 
theoretical models, and both theoretical models and experiments predict the same 
trend in aspect ratio: the aspect ratio increases with laser power. 

Masud, Kamotani, and Ostrach [Ref. 19] experimentally investigated high 
Prandtl number flow in a half-zone configuration. Specifically, they studied the ef¬ 
fects of buoyancy and column shape on the onset conditions of oscillations. There 
is a minimum critical Marangoni number, below which no oscillations appear. They 
speculated that both convection and a dimensionless surface deformation parameter 
must be taken into account together to explain the oscillation mechanism. Kamotani 
and Ostrach [Ref. 20] then conducted a theoretical analysis of the half-zone con¬ 
figuration. They concluded that free surface deformation plays an important role 
in oscillatory thermocapillary flow in high Prandtl number fluids. The deformation 
induces a chang e in the surface flow and alters the driving force in the hot corner. 
This then triggers oscillation cycles in which the flow at the surface alternates be¬ 
tween fast and slow. Finally, they derive a surface deformation parameter by scaling 
analysis, alluded to in their earlier paper [Ref. 19]. The scaling for the core velocity 
found by Cowley and Davis was confirmed for several high Prandtl number flows in 
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the Kamotani and Ostrach paper. Kuhlmann [Ref. 3] points out that this behavior 
suggests that for high Prandtl numbers and Reynolds numbers that axe not too large, 
the hot corner can determine the scalings of the surface temperature and the core 
velocity. 

As stated earlier, the Marangoni effect is caused by a jump in shear stress 
across the interface which balances the surface-tension gradient along the interface. 
Scriven and Sternling [Ref. 13] give a short but interesting historical review of 
this phenomenon, including causes and resistance to surface movements. Batischev, 
Kuznetsov, and Pukhnachov [Ref. 21] investigated the Marangoni effect, dividing 
the cases involving the thermocapillary effect into classes for which the motion shows 
similarity, but they mainly considered the case where the Reynolds number is high 
and the Prandtl number is of order unity. They define the Marangoni boundary layer 
as a surface layer where intensive convection quickly decays in depth. The fluid in 
the Marangoni boundary layer differs in a qualitative sense from that in the Prandtl 
classic boundary layer, which is initiated by an external fluid flow due to the no-slip 
condition at a rigid wall. The Marangoni boundary layer is initiated by forces acting 
on a free surface and subsequently initiates the fluid flow inside the entire volume [Ref. 
21]. Pukhnachov [Ref. 22] systemized the similarity solutions of the Navier-Stokes 
equations in another paper, distinguishing and justifying boundary-layer asymptotics 
when the viscosity goes to zero, describing the construction of nonstationary bound¬ 
ary layers on a liquid free surface and providing descriptions of asymptotic solutions 
in several cases. 

Napolitano [Ref. 23] also studied Marangoni boundary layers between two 
immiscible fluids in the plane, restricting analysis to the case with a Prandtl number 
of order unity. Three necessary conditions for the existence of a Marangoni boundary 
layer were given. First, the region must be sufficiently far away from solid boundaries; 
second, surface driving forces must be of the same order of magnitude as viscous forces; 
and third, the thickness of the region must be much smaller than the thickness of the 
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interface between the two immiscible fluids. The influence of the flow field on the 
dynamic interface’s shape was also calculated in terms of characteristic length ratios. 
Napolitano [Ref. 23] and Ostrach [Ref. 24] separately derived the scaling law of low 
Prandtl number inertial flow in a free surface layer. 

Strani, Piva, and Graziani [Ref. 25] also studied thermocapillary convection 
in a rect angular cavity. They found an asymptotic solution in the limiting case where 
the aspect ratio goes to zero. For increasing values of the aspect ratio, the motion 
is confined in a region near the free surface. As with other investigators, Strani, 
Piva, and Graziani concluded that any surface deformation seems to have a negligible 
influence on the qualitative aspects of the flow field. The solution method employed, 
although independently derived, turns out to be a special case of the method proposed 
by Sen and Davis [Ref. 11]. Regarding the cold comer, they concluded that the larger 
surface temperature gradient near the cold wall develops a local increase in the driving 
force. This leads to a complex stagnation flow field, with significant acceleration near 
the wall, positive until the maximum velocity is reached, then negative to meet the 
boundary condition of zero velocity. A boundary-layer-type flow grows downward 
along the cold wall, starting from the stagnation point at the surface. 

Canright [Ref. 1] introduced study of the dynamics restricted solely to the cold 
corner region, where the flow along the free surface toward the cold wall compresses 
the thermal gradient, thereby enhancing the flow in a sort of positive feedback. This 
is known as the cold corner problem. This feedback results in small local length 
scales and high velocities near the corner. Numerical results are consistent with 
scaling analysis. In particular, for frilly convective flow, the comer behavior is lo¬ 
cally deter min ed. Scaling limi ts for the Reynolds and Marangoni numbers are given 
for four different regimes: conductive-viscous, convective-viscous, conductive-inertial, 
and convective-inertial. Canright concluded that increasing the global Marangoni 
number decreases the local length scale to give an effective local Marangoni number 
of unity [Ref. 1]. In contrast to the hot corner problem posed by Cowley and Davis, 
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Canright found that no thermal boundary layer forms when the Prandtl number is 
large. This is due to the surface forcing being limited to a relatively concentrated 
region in the cold corner, while for the hot corner, the forcing is distributed over a 
broad region, since the thermal variations in the horizontal direction are extended by 
convection. Canright also validated his scales by comparing his results with those of 
Zebib, Homsy, and Meiburg [Ref. 17], who showed that the temperature gradient on 
the free surface can vary considerably within a small distance from either of the hot 
and cold walls. Canright’s scaling of the boundary layers were also published in a 
treatment of the cold corner by Kuhlmann [Ref. 3]. 

Huber [Ref. 26] extended the study of the cold corner problem, examining 
the problem numerically for different low Marangoni numbers using a Green’s func¬ 
tion flow representation for the viscous case, in the limit as the Reynolds number 
approaches zero. One advantage of the Green’s function is that the flow can still be 
represented over the entire quarter plane, so that there is no artificial recirculation 
due to imposed artificial boundaries. The Green’s function approach does not require 
flow boundary conditions at the computational domain’s boundaries. The flow is 
assumed isothermal across the boundary and is recirculating, decaying with distance. 

Canright’s work in 1994 concentrated on the Marangoni number dependence 
of the cold corner problem. A representative velocity vector field (using Canright’s 
numerical data) is shown in Figure 1. This work reconsiders the problem, re finin g 
the Prandtl number dependence of the variables. All the graphs shown in [Ref. 1] 
are correct, but we will now show that the theoretical scales in the cold corner for 
the convective inertial case are different than those previously published, in terms 
of the Prandtl number dependence. Canright’s numerical results reflect a boundary- 
layer structure that we assume throughout this work and use as a guiding principle. 
The task is to solve for the flow and/or heat in each region and match it to that in 
neighboring regions in order to predict future flow conditions. The scalings derived 
herein are different than those published by Canright; nonetheless, they agree better 
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Distance From Free Surfa< 



Figure 1. Velocity Vector Field 


with his numerical data. 
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III. PROBLEM STATEMENT 


Most of the fundamental ideas of science are essentially simple, and may, as 

A RULE, BE EXPRESSED IN LANGUAGE COMPREHENSIBLE TO EVERYONE. —Albert Einstein 


A. MATHEMATICAL FORMULATION OF THE PROB¬ 
LEM 

A pool of incompressible Newtonian fluid is bounded on the left by a vertical 
solid wall, maintained at a cold temperature T c , while the undisturbed fluid far from 
the cold corner is at the hot ambient temperature Th of the core flow (see Figure 2 
below). 

AIR 



Figure 2. Problem Formulation 


Above the horizontal free surface of the liquid is an inviscid, nonconducting 
gas. Surface tension is assumed strong enough to keep the free surface flat (small 
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Capillary number), but with surface tension variations due to a linear dependence on 
temperature. The resulting flow is assumed to be two-dimensional and steady. 

The equations governing the thermocapillary convection in the cold corner are 
the conservation of mass, the conservation of momentum, and the conservation of 
energy equations: 

V • u = 0 (III.l) 

pu-Vu = -Vp + pV 2 u (III.2) 

pcpU- VT = kV 2 T (III.3) 

with the boundary conditions that at y = 0, T y = 0, v — 0, and pu y = ^T x . At 
x = 0, T = T c , and u = v = 0. As (x, y) —► oo, T —> Th, and u, v —> 0. 

In these equations, u is the velocity vector with components u and v in the 
x (horizontally rightward) and y (vertically downward) directions, p is the pressure, 
T is the temperature, p is the density, p, is the viscosity, Cp is the specific heat, k 
is the thermal conductivity, and the surface tension is assumed to be of the form 
<7 = er c — 7 (T — T c ), with 7 a positive constant. The boundary conditions specify 
that the cold wall is isothermal with fluid obeying the no-slip condition, and the flat 
free surface is thermally insulated, with thermocapillary forcing. 

To nondimensionalize the equations, a heat flux scale of Q and a temperature 
scale (relative to the cold temperature) of AT = Th~T c axe introduced. Thermal con¬ 
duction gives the length scale of d= Q/k AT, the thermocapillary coupling gives the 
velocity scale u = 7 A T / p, and the viscous pressure scale is p = pu/d = k-y (AT) 2 /Q. 
In summary, the nondimensional scales are: 


X = 

(Q/k AT) T 

(in.4) 

u = 

{l AT/p) u' 

(III.5) 

p = 

(k 1 (ATf/Q)p > 

(III.6) 

T = 

(T-T h )/(T C -T h ) 

(III.7) 
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which implies T — T^ — (AT) T 1 , where the variables with primes are dimensionless. 
Begin with the conservation of energy equation, which can be rewritten as 

u- VT = kV 2 T (III.8) 

where k is the thermal diffusivity. Substituting the dimensionless scales into this form 
of the energy equation yields 

^[u'.VT'] = V 2 r (III. 9 ) 

Dropping the primes and defining a dimensionless parameter M, the Marangoni num¬ 
ber, asM = ^|, the nondimensional conservation of energy equation becomes 

M(u-VT) = V 2 T. (III. 10) 

Next, the conservation of momentum equation may be rewritten as: 

u ■ Vu = —^ Vp + v V 2 u (III. 11) 

where v is the kinematic viscosity. Substituting the dimensionless scales into the 
momentum equation yields 

p • W] = -Vp' + V 2 u'. (III. 12 ) 

Dropping the primes and defining a dimensionless parameter R, the Reynolds number, 
as R = the nondimensional conservation of momentum equation becomes 

R (u • Vu) = -Vp + V 2 u. (III. 13) 

Situations in which the Reynolds number is small are called slow viscous flows, because 
the viscous forces arising from shearing motions of the fluid predominate over inertial 
forces associated with the acceleration or deceleration of the fluid particles. The 
thermocapillary boundary condition is also nondimensionalized by substituting (III.4) 
and (III.5) into pu y = 7 T x , obtaining (dimensionless) u y = T x . In addition, the mass 
conservation equation, V • u = 0, remains unchanged after nondimensionalization. 
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The ratio of the Marangoni number and the Reynolds number gives the nondi- 
mensional Prandtl number 

P = v/k = M/R. (III. 14) 

We want to analyze the system for large values of the Marangoni number, M, which 
measures the importance of thermal convection relative to thermal diffusion. For large 
enough M, convection becomes important and the strong surface flow towards the 
cold wall compresses the thermal gradient along the surface, which in turn strengthens 
the local driving force for the flow. Also, since R = M/P , the effect of making R 1 
is to make the inertia force much greater than the viscous force, so that inertia and 
pressure forces are dominant, except along the rigid wall and free surface, where 
viscous boundary layers are formed. 

Near the surface, vorticity must be taken into account. Vorticity is the curl 
of the velocity, Q = V x u, which (in two-dimensional flow) reduces to the one 
component along the axis normal to the x-y plane, yielding u = v x — u y . Thus, 
eli min ating pressure by taking the curl of the momentum equation, the following 
equation is obtained: 

u • Vu) = v V 2 u> (III. 15) 

This equation is referred to as the vorticity transport, or vorticity transfer, 
equation. It states that the substantive variation of vorticity, which consists of the 
local and convective terms, is equal to the rate of diffusion of vorticity through friction. 
In our non-dimensionalized form, the vorticity equation becomes R (u ■ Vo;) = V 2 u. 
At the surface, the thermocapillary condition gives uj = T x , which comes from u = 
v x — u y , and the fact that v x = 0 at the surface. 

B. ASSUMPTIONS AND SCALING ANALYSIS 

A few assumptions concerning the scaling analysis are made. First, a high 
Reynolds number, R, implies that inertial forces are dominant, outside the viscous 
boundary layers. Second, a high Marangoni Number, M, implies that thermal con- 


18 



vection is important. Third, a low Prandtl Number, P, implies that there exists a 
low ratio of viscous to thermal diffusion and the region is inertial. 


AIR Ty = o 



Figure 3. Scaling of the Domain 


Recall that we are assuming the structure based on Canright’s numerical data. 
Figure 3 shows a schematic of the corner with viscous boundary layers along the free 
surface and rigid wall. The depth and width of the pool are large compared to all 
local length scales, so the pool appears semi-infinite both horizontally and vertically. 

To distinguish the various notations, a hat Q will denote the original nondi- 
mensionalized variables. Capital letters will denote the core scaling, and lowercase 
letters will denote boundary-layer variables that do not have the same scaling as the 
core variables. 

Let the core velocity scale be denoted by U. Let the characteristic velocity 
scale along the free surface be u. We assume that surface thermal variations will 
be confined to a region with a horizontal characteristic length scale of l. We also 
assume that vorticity will be confined to thin viscous boundary layers. The flow in 
the core region is dominated by inertia. (This “defines” the core region as the area 


19 





with negligible vorticity). The solid boundary acts as a source of vorticity, which is 
then diffused away by viscosity and converted downstream with the fluid. 

Along the surface, the boundary layer results from thermocapillary stresses 
and will have thickness 8. Along the wall, the boundary layer results from the no-slip 
condition and will have thickness A. Thus, 8 is the vertical length scale of velocity 
shear along the free surface, and A is the horizontal length scale of velocity shear 
along the rigid wall. 

In case of large M, convection is important and the rapid surface flow into 
the cold corner compresses the thermal gradient along the surface, reducing l. The 
length l will be reduced to the point that thermal diffusion away from the rigid wall 
will balance the strong convection toward the wall. 

In case of large R, inertia is dominant and there will be two viscous boundary 
layers of thickness 8 and A (~ 8) along the surface and rigid wall, respectively. The 
effect of increasing the Reynolds number to values large compared with unity is to 
confine the vorticity diffused from the rigid wall to a layer of relatively small thickness. 
As R increases, 8, the thickness of the boundary layer, decreases. Convection contains 
the vorticity inside the boundary layer. 

Surface tension variations are due to an assumed linear dependence on the 
temperature. For the boundary layer along the free surface, the thermocapillary 
stress condition at the surface gives u y ~ T x . This provides the relation | ~ or, 
solving for the surface velocity in the horizontal direction, u ~ f. The continuity 
condition V • u = 0 gives the relation u x ~ v y . Therefore, y ~ or v ~ j, which 
gives v ~ (j) 2 as the vertical component of the velocity near the surface. 

At the boundary-layer edge/core region interface, normal velocities must match, 
so v ~ V, and continuity requires that V y ~ U x . Thus V y ~ ^ ~ U x ~ j, which 
gives the relation ^ ~ f • 

Next, a check of the dominant balance of terms is required. In the vorticity 
equation, R (jp +Ignore ^ as too small, since l 3> 8. This implies 
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that R J| ~ jfi, or that R~- ~ 1 . Therefore, this gives an expression for u in terms 
of the Reynolds number, w,v ^>or«~^. 

Similarly, in the energy equation applied to the core region, MUT X ~ T xx (this 
dominant balance will be described in detail later; for now, assume T y and T yy scale 
comparably to T x and T xx in the core flow region). Thus, MUT X ~ Mj ~ T xx ~ 
so that ~ p, which yields another relationship for the surface velocity, u ~ 

We assume that the wall boundary layer is passive, fed by the mass flux from 
the surface layer (due to conservation of mass). We can derive its scales by using the 
mass flux, continuity and vorticity equations. Conservation of mass implies that the 
velocity out of the corner is the same as that into the corner (v ~ u). The moment um 


equation gives R(Uv x + vvy) — v xx + vyy, and continuity gives U x ~ vy. This means 
a ~ Z> 01 v ~ Substituting, ^ ~ » xp, leaving 

A ~ M~ l PU~ l and L ~ 


To obtain scales for the lengths and velocity, use the three relationships for 
the free surface horizontal velocity scaling: Solving for the 

horizontal boundary layer length l, l ~ 6 2 M ~ p, which implies that <5 ~ j±p. Thus, 
l ~ 7 ^ 2 , u ~ P, and U ~ P 2 . Further, A ~ ~ 6 and L ~ ~ /. 

To summarize, using the boundary condition at the free surface gives a starting 
point for the scalings and dominant balance. It is important to obtain relations for 
the thermal length scale l, the viscous thickness 6 , the surface velocity u, and the 
core velocity U scales in the cold corner, in terms of the dimensionless parameters M 
and P. These scalings are 


6 ~ M - 1 P _1 , 
l ~ M~ l P ~ 2 , 
u ~ P, 

U ~ P 2 . 
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These scalings are consistent with the numerical results obtained by Canright in his 
1994 work [Ref. 1]. 
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IV. THE COLD CORNER REGIONS 


It is not knowledge, but the act of learning, not possession, but the act of 

GETTING THERE, WHICH GRANTS THE GREATEST ENJOYMENT. —Carl Friedrich Gauss 


A. OVERVIEW 

The flow into the corner is assumed to be convective and inertial, due to the 
large Marangoni number and small Prandtl number (M » 1 ,P <C 1), governed by 
the system of equations 

R (Z-Vu) = V 2 £ (IV. 1) 

V-fi = 0 (IV.2) 

M (Z-VT) = V 2 f (IV.3) 

The cold corner is divided into four regions: an outer core region away from 

the surface and cold wall where the flow is inviscid and relatively simple; an inner 
viscous boundary layer along the liquid free surface; an inner viscous boundary layer 
along the wall/liquid interface; and a thermal layer overlapping the viscous layers and 
the core flow. A representative velocity vector field (using Canright’s numerical data) 
was shown in Figure 1 in Chapter II. The task is to solve for the flow and/or heat 
in each region and match it to that in neighboring regions in order to predict future 
flow conditions. The scaling factors are 6 ~ A ~ M~ l P~ l , l rv Jj M~ l P ~ 2 , u ~ P, 
and U ~ P 2 . In the core flow region, let the horizontal and vertical components of 
the velocity be designated by capital letters. The governing equations in each region 
will be examined, using the appropriate nondimensionalized scalings. 
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B. VISCOUS SURFACE BOUNDARY LAYER 

A viscous boundary layer forms along the surface, with thickness 6 and length 
l. The scalings used in this region are u = Pu, v = P 2 V, x = M~ l P~ 2 X, y = 
M~ l P~ l y, T = T + • • •, and p = M P z p. The pressure scales with the dominant 
term in the core flow. The momentum equation in the x-direction can be written as 

UUX + Vu y = ~P 2 PX + P 2 + Uyy, (IV-4) 

where to leading order we can neglect the terms of order P 2 . The boundary conditions 
at y = 0 give u y = Tx and V = 0. The continuity equation becomes 

u x + Vy = 0. (IV.5) 

Prom the momentum equation in the y-direction, the pressure change in the y- 
direction is negligible through the boundary layer, indicating that the pressure is 
effectively constant throughout the layer for each X value. Mathematically, px p y - 
It is assumed then that the pressrue distribution in the boundary layer is imposed by 
the core flow pressure gradient just outside the boundary layer. In Ludwig Prandtl’s 
words, “the pressure distribution will be impressed on the transition layer by the free 
fluid” [Ref. 27]. Also, the velocity in the free surface boundary layer is continuous 
with the core velocity U. 

The energy equation in the surface viscous boundary layer is derived using a 
boundary-layer expansion for T and will be given later. 

C. VISCOUS WALL BOUNDARY LAYER 

A viscous boundary layer forms along the wall, with thickness A and length 
L. This viscous boundary layer is different from the surface viscous boundary layer 
in that there is a no-slip condition at the wall. The scalings used in this region are 
u = P 2 U,v = Pv,x — M~ l P~ l x, y = M~ l P~ 2 Y, T = T + • • •, and p = M P 3 p. 
The momentum equation in the y-direction can be written as 

U V x + V Vy — — P 2 py + V xx + P 2 Vyy. (IV.6) 
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where again terms of order P 2 can be neglected to leading order. The no-slip condition 
at the rigid wall gives u = v = 0atx = 0. Similarly, the continuity equation becomes 

U x + v Y = 0. (IV.7) 

The energy equation for the viscous boundary layer on the wall is 

UT x + vT y = P~ 1 T xx + P T YY (IV.8) 

with the boundary condition at x = 0 of 

T = -1 (IV.9) 

The PTyy term is neglected as small. The T xx term must therefore be small enough 
to balance the equation. Thus, there is no need to solve the energy equation in this 
region, as the temperature is basically constant across the wall’s viscous layer. Also, 
due to the conditions of high Marangoni number coupled with low Prandtl number, 
the larger thermal layer encompasses the viscous boundary layer at the wall. 

D. CORE FLOW REGION 
1. Flow Equations 

The core flow is assumed to be irrotational, which means that there is no 
vorticity, and also is incompressible, which means that density is not affected by 
changes in the pressure. When the velocity U is solenoidal, we may introduce a 
stream function ip which allows us to replace the horizontal and vertical components 
of velocity by a single function. Define ip by u = and v = —1|. These defin¬ 
itions satisfy the continuity equation. When the flow is also irrotational, then the 
expression of zero vorticity gives Laplace’s equation, V 2 ip = 0. Although the equa¬ 
tion of motion is nonlinear in U, the velocity distribution is completely determined 
by a linear equation derived from the restrictive condition of irrotationality and the 
mass-conservation equation. The nonlinear equation of motion is needed only for the 
calculation of the pressure after the velocity distribution has been determined. The 
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boundary conditions for the core flow must match the normal velocities in each of 
the boundary layers. Scale the velocity in this region by P 2 . Scale the horizontal and 
vertical lengths by M~ l P~ 2 . Scale the pressure by MP 3 . The momentum equation 
in the ^-direction is 


U Ux + VUy = —px + P ( Uxx + Uyy)- (IV. 10) 

The momentum equation in the y-direction is 


U V x + V V y = -Py + P(V X x + V YY ). (IV. 11) 


With M 1 , and P « 1 , these momentum equations reduce to U Ux + V Uy = ~Px 
in the ^-direction and UVx + VVy = ~Py in the y-direction. The momentum 
equations could also be written in terms of the stream function, ip. For example, the 
momentum equation in the ^-direction becomes 

1p Y IpXY ~ '’Px 1 pYY = — PX + P(' l PxXY + IpYYY ) (IV. 12) 


Laplace’s equation can be written as 


V 2 ip 


d 2 ip d 2 ip 

dX* + dY* 


(IV. 13) 


In a two-dimensional flow, the lines of constant ip are known as the streamlines, and 
the difference between the numerical values of two streamlines is equal to the flow rate 
between the streamlines. An equivalent description of two-dimensional incompressible 
irrotational flow is in terms of a velocity potential (p , such that u = V</>. The stream 
function ip is harmonic, as is 0, and together they satisfy the Cauchy-Riemann equa¬ 
tions as conjugate harmonic functions, which means that lines of constant ip and (p 
are orthogonal. The fact that ip and (p are harmonic and satisfy the Cauchy-Riemann 
equations is necessary and sufficient for the definition of a complex potential F, de¬ 
fined by F = (p + iip, where z = x + iy and the stream function and velocity potential 
are both functions of x and y. F is analytic, and note that F'(z) = (p x (x, y ) +iip x (x, y), 
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or F'(z) = <fi x (x, y) - i<p y (x,y). The velocity thus becomes u = F'(z), and the speed, 
or magnitude of the velocity, is found by \u\ = \F'(z)\. 

Since the equation V 2 ip = 0 is linear, we may superpose solutions for different 
flows and add directly the values of ip at every point in the plane to obtain the new 
values of ip, which represent physically the direct superposition of the various flows. 

2. Dominant Balance 

Naive APPLICATION OF THE METHOD OF DOMINANT BALANCE IS BOOTLESS. — Carl M. 

Bender and Steven A. Orszag 

How does the temperature T change across the boundary layer? Let U be the 
core velocity and T be the core temperature, and u and t are the additional values 
of the velocity and temperature inside the surface viscous boundary layer, so u and 
t are significant only inside the layer. The horizontal derivative of the temperature, 
T x , in the core flow region is similar to T x at the surface. 

Let T = T(X,Y) + rt(X,y), where r is an unknown scale. Further, let 
u = P 2 U(X, Y) + Pu(X,y), using the known scales: x = y = ^ = JL. 
To determine the size of r, substitute the expressions for T and u into the energy 
equation. This yields 

M [( P 2 U + Pu)(T + Tt)x + P 2 V(T + rt) t ] = (T + rt )« + (T + rt)yy (IV. 14) 

or, in expanded form, 

M 2 P a UT x + M 2 P 3 uT x + M 2 P i Urt x + M 2 P z urt x + M 2 P*VT Y + M 2 P 3 Vrt y 

= M 2 P 4 T xx + M 2 P 4 rt xx + M 2 P a T yy + M 2 P\t yy 

Divide this equation by M 2 P 4 and subtract the core flow region equation, UT X + 
VT y = T xx + T yy , which leaves 

pUT x + rUt x + —rut x + -prVty = rt xx + rt yy . (IV.15) 
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What is the do min ant, balance? That is, how must the scale r be chosen to balance 
the additional heat convection of the first term in (IV. 15)? 

Suppose p uT x ~ rUt X - Then r ~ V, which implies the next term -prut x is 
0(p 2 )- Since that is much larger, the assumed balance cannot be dominant. Thus, 
the assumption is not consistent. 

Suppose puT x ~ prutx■ Then r is 0(1), and the term p^rtyy is larger, which 
means the assumed balance is not dominant and the assumption is not consistent. 

Suppose puTx ~ prVt y . This is equivalent to the previous case, so it is not 
consistent. 

Suppose -puTx ~ Ttxx- Again, this would mean that r ~ p, as in the first 
case, which is not consistent. 

Suppose puT x ~ -j&Ttyy. Then r ~ P and all the other terms are small in 
comparison to the dominant terms. This is the only consistent case. Hence r — P. 

Thus, the change in temperature T across the boundary layer is Pt, i.e. O(P), 
so it is negligible. This implies that the temperature gradient Tx at the top of the 
core is equivalent to that right on the surface, which drives the boundary-layer flow. 
But the vertical temperature gradient Tg = MP 2 [TV + t y ], so the surface boundary 
condition Tg = 0 implies TV = —t y at y = 0. The heat equation for the surface 
boundary layer is 

UT X = tyy. (IV. 16) 

In order to determine how surface layer convection modifies heat flux into the core, 
use the boundary condition Tg = 0 and integrate the left-hand-side of the surface 
heat equation across the surface boundary layer. As y —* 00 , t y —> 0 and the heat 
flux condition in the core becomes 

T Y (X,0)= f°° uT x dy. (IV. 17) 

Jo 
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Since / 0 °° u dy is just the mass flux into the corner (or the expression for the stream 
function, xp), the heat flux condition becomes 


T Y (X,0)=T X 7p(X,y->oc). 


(IV. 18) 
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V. APPROXIMATE METHODS 


An idea which can be used once is a trick. If it can be used more than once 

IT BECOMES A METHOD. —George Polya and Gabor Szego 


The goal of this chapter is to match the velocity in the surface viscous boundary 
layer to that in the wall viscous boundary layer and to that in the core flow region. 
The temperature can then be found using the core velocity. For the viscous boundary 
layers, several approximate methods exist based on integrated forms of the boundary- 
layer equations. These are integral methods which do not attempt to satisfy the 
governing equations of the boundary layer for each streamline. Instead, the equations 
are satisfied on average over the thickness of the boundary layer. First, start with 
the viscous boundary layer along the free surface, since the surface thermal gradient 
drives the flow into the comer. 

A. VISCOUS SURFACE BOUNDARY LAYER 
1. Background 

The first approximate method considered is that developed by Pohlhausen in 
1921 [Ref. 28]. It is based on integrating the momentum equation across a boundary 
layer. The integrated momentum equation gives an ordinary differential equation 
for the thickness of the boundary layer, S(x), provided that a suitable form for the 
velocity profile is assumed. Rosenhead [Ref. 29] and Schlichting [Ref. 30] offer two 
excellent overviews of this method, as well as descriptions of updated methods in more 
modern forms. The key to the approximate methods is that they allow calculation 
of the momentum thickness, the displacement thickness, and the shearing stress at 
a wall. These will all be defined later. Most of the approximate methods for two- 
dimensional flows in the literature deal with flow past a rigid wall. No such method 
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for flow along a free surface (where the no-slip condition of the rigid wall does not 
apply) has yet been discovered in the literature. 

What is the critical boundary-layer thickness, 8 (x)? Pohlhausen’s method de¬ 
termines the boundary-layer thickness by first writing the momentum equation in 
terms of the shear stress, the displacement thickness, and the momentum thickness 
of the boundary layer. The displacement thickness 81 has a simple physical interpre¬ 
tation. U 6 \ is the decrease of the mass flux across a normal to the surface, due to the 
boundary layer. The stre amlin es of the outer flow are thus displaced away from the 
surface through a distance 6 \ [Ref. 29]. It is the distance that a wall would have to 
be displaced outward into the free stream in order not to change the flow field if the 
fluid were completely inviscid and there were no boundary layer. The fluid becomes 
slowed down by viscosity near the surface and is thus forced outward, so that the 
effective surface presented to the oncoming stream is increased (or thickened) by the 
boundary layer’s displacement thickness [Ref. 32]. The momentum thickness has the 
same physical significance except that it is based on momentum flux instead of mass 
flux. 

Pohlhausen’s approach introduces a dimensionless quantity A, a shape factor, 
which can be interpreted physically as the ratio of pressure forces to viscous forces. 
According to Schlichting [Ref. 30], in order to obtain a quantity with real physical 
significance, the thickness 8 must be replaced by a length which itself has physical 
significance, such as the momentum thickness. The fundamental idea of approximate 
methods based on the momentum thickness is to assume that the dependence of 
the solution, u(x, y), on y can be expressed by some known expression of y, in which 
coefficients appear, which are unknown functions of x. As the value of y gets large (the 
distance from the free surface increases), the boundary-layer velocity u approaches 
the core velocity U outside the boundary layer, and its derivatives tend to vanish. 
Generally, only as many boundary conditions for y = 0 are taken into account as is 
necessary to determine the unknown coefficients in the velocity profile. This will be 
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done later in this section. 

In view of this, a more accurate method for the present situation is that of 
Timman [Ref. 31], whose method is a variant of Pohlhausen’s. In 1949, T im man 
developed a family of profiles to model the velocity near a no-slip surface more ad¬ 
equately than Pohlhausen. Timman’s method seeks to find an exponential form for 
the velocity profile u(x, y) which satisfies the momentum equation and some of the 
boundary conditions, in a case where the velocity profile exponentially decays through 
the boundary layer. As Rosenhead states, “It is hoped that this form will approximate 
to the exact profile, which satisfies all of the conditions” [Ref. 29]. The complimen¬ 
tary error function which is introduced by Timman produces the desired asymptotic 
behavior; namely, the velocity exponentially decays as the distance from the surface 
increases. We will attempt to use a modification of Timman’s method in the case of 
the velocity flow along the free surface. 

2. Timman’s Method 

The velocity profile form assumed is 

§ = /M. ” = ski < V1 > 

where <5(X) is the effective total thickness of the boundary layer and U(X) is the 
velocity of the external flow. The function / may also depend on X through some of 
the coefficients which are chosen to be part of the profile in order to satisfy surface 
boundary conditions. Timman’s profile is of the form 

jj = f(r)) = l-J e- 7 ? (a + erf-\ -) drj — e _7?2 (6 + drf H-), (V.2) 

where a, b, c, and d, etc., are coefficients which will turn out to be functions of X. 
As y —» oo, it follows that / —► 1 to match the external flow. Therefore, only the 
boundary conditions at the no-slip surface rj = 0 are considered. Only the same 
number of boundary conditions as unknown coefficients are used. 

In Timman’s profiles, the velocity difference f(rj) — 1 = ^ — 1 exponentially 
decays through the boundary layer and satisfies u = 0 at the surface. In the core 
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region, 77 —> 00 , f(rj) —> 1 , and u —> U, or the velocity in the boundary layer matches 
the core velocity. 

In our problem, the surface is a free surface, and the no-slip condition of Tim- 
man’s method does not apply, since 0 at the surface. The velocity in our modifi¬ 
cation must be adjusted accordingly to match the free surface boundary conditions. 
The velocity u stills decays exponentially through the boundary layer, but where Tim- 
man’s approach assumes ^ = 0 at the surface and ^ = 1 at the core/boundary-layer 
edge (where the boundary layer meets the core flow), our problem must be modeled 
to have u as a maximum at the surface and u—U —* 0 at the core/boundary-layer 
edge, since from Canright’s numerical data [Ref. 1 ] and the scaling analysis, the core 
velocity is negligible. 

Consider again the momentum equation in the free surface boundary layer: 

UU X + Vu y — Uyy (V-3) 

with boundary conditions of V = 0 and u y = Tx when y — 0. Combine u times the 
continuity equation with this momentum equation: 

U ( U X + Vy = 0) © (UU X + Vu y = Uyy ) 

which gives 

UU X + UVy + UU X + Vu y = Uyy 

This can be rewritten as 

Nr + (uV)y = ( Uy)y . (V-4) 

Integrate this expression with respect to y and simplify: 

rOO pOO pOO 

/ (uu) x dy+ / (uV)ydy = / (u y ) y dy (V.5) 

Jo Jo Jo 

poc 100 100 

2 Jo uu * d y + uV \0 = %| 0 = -Wy=0 ( v - 6 ) 

Since u = V = 0asy— >00 and V = 0 at y = 0, the expression for u y along the 
surface becomes 

iV“y) y =o = 2 [ uu x dy- (V.7) 

JO 


34 



To begin the approximation, use a value for Tx = u y , the free surface boundary 
condition. Along the surface, Tx decreases sharply in magnitude as the distance from 
the wall increases, decaying along the length of the boundary layer. This can be seen 
in Figure 4 , which shows the gradient of the surface temperature in the r-direction, 
using Canright’s numerical temperature data. It also shows the temperature profile 
(dotted line) along the surface. 

Can the revised momentum equation be put into a form that has a displace¬ 
ment and momentum thickness? Expand the upper limit of integration in the mo¬ 
mentum equation to infinity, to make it consistent with the range of u(X, y). Define 
62 as the momentum thickness of the boundary layer. Let 

62 = J u 2 dy (V.8) 

Since 77 = |, where S(X) is the surface boundary-layer thickness, y = Sr), and dy = 
6 dr). 

If }(rj) is defined by f(rj) = u, then 

d 2 = S f f 2 {rj) dr). (V.9) 

J 0 

Thus, the skin friction, or shear stress, can be written in terms of the momentum 
thickness, yielding a first-order differential equation r 0 = — ^ (<5 2 ), or 

JL(6 2 ) + t 0 = 0 (V.10) 

3. Modification to Timman’s Method 

Consider a family of velocity profiles where all of the coefficients are functions 
of X, as shown in (V.2). If no terms of higher order than the c(X) and d(X) terms 
are included, then the modification to Timman’s profile yields 

r oo 

u(X, y) = f(r)) = J e~ t2 (a(X) + c(X)t 2 ) dt - e^ 2 (b(X) + d(X)r) 2 ) , (V.ll) 
where 77 = g(xj‘ Wh en y = 0, the surface velocity becomes 

u(X, 0 ) = ^a(X) + ^c(X) - b(X). (V.12) 
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Tx 



Figure 4 . Temperature Gradient along the Free Surface 


Due to the complications of this model, we simplify it even further (see Ap¬ 
pendix A for our attempt at deriving a full Timman-type model): let b(X) = c(X) = 
d(X) = 0. Then 

u = f(v) = f°° a(X)e~ t2 dt, (V. 13 ) 

Jrj 

or 

f(v) = [1 - erf(77)] = erfc(77), (V. 14 ) 

where erfc(77) is the complementary error function of 77. A typical characteristic of 
the complementary error function of 77 is that as 77 increases, erfc(77) quickly decays 
to 0. 

Solving the ordinary differential equation (V. 10 ), ^ (S 2 ) + To = 0 , for the 
momentum thickness yields 

S 2 = - f X To dX' + C\ (V. 15 ) 

Jo 
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However, the shear stress tq = u y = T x from the given surface boundary condition, 
so 


62 


[ X T X 'dX' + c u 
Jo 


or 62 = — T + ki. The boundary condition that must be satisfied is that as X —» 00 , 
82 (X) —> 0, u —> 0, T —0, and tq —» 0. This implies that k\ = 0. Thus, 


62 = -T. (V.16) 

(Since the temperature is less than zero, 5 2 is a positive quantity). The momentum 
thickness has a value equal to the magnitude of the temperature. This implies two 
things. First, the choice for the velocity profile does not explicitly affect the solution 
of the momentum thickness integral. Second, since the temperature is of order 1, 
the momentum thickness is of order 1. The question that remains is: what is the 
boundary-layer thickness 81 

Since f(rf) = |a(X) y / 7rerfc(r;), substitute this expression into the equation 
for the momentum thickness (V.9) and obtain 


82 = -T 


(V.17) 


- 


[erfc(77)] 2 dr] 


(V.18) 


Now determine 8 (X) in terms of the temperature. First, the coefficient a(x) 
needs to be determined. Recall that f(rj) = u, so that 


u 


y ~ 


df 

dy 


1 

8 (X) 




where the prime denotes differentiation with respect to rj. Thus, f(rj) = 8 (X) u y = 
8 (X)T x along the surface. 

Differentiating the velocity profile yields 


/'(*?) = ~a(X) e~< 


(V.19) 
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Along the surface, y = r\ = 0, so that /'(0) = —a(X), which implies a(X) = 
- 6 (X) u y = - 6 (X) T x . Thus 

r r2 nr rOO 

62 = S 3 (X) / [erfc(7/)] 2 drj = -T (V.20) 

4 Jo 

This allows the determination of the boundary-layer thickness S(X) to be 

rp 

HX)= -rTZ ■ (V.21) 

. 1 X £ . 

where 

£= i L f erfc ^ 2 dr] = ( : _ ^r) • ^ V - 22 ^ 

Numerically, e « 0.2596 (see Appendix B for derivation). So, 6 (X) is positive and of 
order 1, which is not unreasonable. 

Assembling the components of the profile, the velocity can be modeled as 

u = — ~~tJ~ 1/3 (-T) 1 / 3 (T x ) 1/3 erfcfa). (V.23) 

Examining each term, (—T) 1//3 decays to 0 as the distance from the wall, X, 
gets large. (Tx) 1 / 3 also decays to 0 as X increases. In addition, erfc(r?) decays 
rapidly to 0 as the distance from the surface, y, increases, which shows that u goes 
to 0 at the edge of the surface boundary layer. As seen in Figure 5, there is excellent 
agreement between the predicted velocity profile along the surface (discrete points) 
and the velocity profile generated from Canright’s numerical data (solid curve). In 
order to evaluate the velocity u, Canright’s surface temperature data was used. Also, 
taking a vertical slice of u through the boundary layer also gives reasonably good 
agreement between the predicted profile and the numerical profile. Typical profiles 
for a y-direction slice of the horizontal component of the velocity, u, are shown in 
Figure 6. As can be seen, the value of u increases exponentially as the depth below 
the surface decreases. In this figure, u increases by more than an order of magnitude 
as it passes through the boundary layer. The predicted velocity is given at distances 
of X = 0.0299 and X = 0.0494 (dotted fines) from the wall, and are compared to 
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Surface Velocity Comparison 



Figure 5. Comparison of Surface Velocity Profiles using Timman’s Method Prediction 
and Canright’s Numerical Data 

velocity profiles taken at the same ^-distances using Canright’s data (solid lines). The 
profiles using Timman’s model decay much faster than the profiles of the n um erical 
data as the distance from the surface increases. This might be due to the lack of 
a second parameter in the model. Also, note that the surface boundary condition 
(u y ~ Tx) matches for each corresponding plot. 

B. VISCOUS WALL BOUNDARY LAYER 
1. Background 

The behavior of the velocity along the vertical rigid wall is similar to that 
outlined in Timman’s paper. The no-slip condition exists at the wall, so the vertical 
component of velocity is zero at the wall, increases to a peak position (where the 
vorticity is zero), and then decreases to match the value of the core velocity at the 
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Figure 6. Typical Profiles of Velocity (u) through the Surface Boundary Layer 

edge of the viscous boundary layer. Two approaches are studied. The first is a method 
by Glauert, which considers the velocity profile through the wall boundary layer as 
similar to that near a plane wall jet. The Glauert approximation matches numerical 
work done by Canright [Ref. 1] in 1994 quite well. A discussion of the vorticity flux 
is also included as an appendix, as the vorticity u can be represented by the gradient 
of the vertical component of the velocity, v x . An interesting conclusion is reached. 
The second method is a Timman-type approximation, but since four parameters 
are needed for accuracy, this method involves solving coupled nonlinear differential 
equations. Since the Timman method worked with the surface boundary layer, its 
use with the wall boundary layer was investigated and is given as an appendix. The 
Glauert approach is more desirable due to its simplicity. 
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2. Glauert’s Method 

Suppose the fluid flow from the free surface into the cold corner turns down the 
wall and acts in the manner of a plane wall jet. There is a similarity solution for the 
velocity profile through the wall boundary layer. Following Glauert [Ref. 36], start 
with the dimensionless momentum equation, assuming the pressure is everywhere 
uniform, 

U v x + v v Y = v xx . (V.24) 

Combine the continuity equation with a definition of a stream function ip, defined 
by U = § and v — — ff. Consider the possibility that there shall be a similarity 
solution of these equations, with the velocity proportional to Y a and the jet thic kness 
proportional to Y b . The two sides of the momentum equation (V.24) will vary with 
Y in the same manner if a + 26 = 1 . For a wall jet, fluid momentum is not conserved, 
due to shear stress from the wall. 

Write the stream function as 


= -Y'-Ofi-n) 


where 


Then 


T) = (1 - b)xY~ b . 


v = (1 - 

and (V.24) becomes 

/'" + ff" + af 2 = 0 (V.25) 

where a = (2b — 1)/(1 — b). The boundary conditions require that /(0) = /'(0) = 0 
and /'(oo) = 0 . Integrating (V.25) between the limits 77 and 00 yields 

roc 

f + ff - (« - 1) / f n dr) = 0. (V.26) 

Jr\ 

Let g(rj) = fjj° f 2 dr), such that 

f + ff-(a-l)g = 0. (V.27) 
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Notice that g > 0. Multiply both sides of the above expression by /', yielding 


/7" + //' 2 -(a-!)»/' = 0. (V.28) 


Integrate again, using integration by parts on the last two terms and thus 

j/' 2 -/» + (<*- 2) jT /'s */ = 0. (V.29) 

If / = 0 when /' = 0, then a = 2, since /“ f'g dr) > 0. This gives a — — \ and b = 

Regarding what Glauert terms the external momentum flux, start again with 
the momentum equation (V.24). Integrate this with respect to x between the limi ts 
x and oo, using the conditions that v —► 0 as x —*■ oo. 



U v x dx + 



vvy dx 



(V.30) 


By continuity, U x = —vy, and from integration by parts, the first term above becomes 

roo roc 

/ U v x dx — —Uv + / vvydx (V.31) 

Jx Jx 


Substituting, 


8 f°° n 

— y v 2 dx- Uv + v x = 0 (V.32) 

Multiply this equation by v and integrate with respect to x between 0 and oo. 

f) roc / roo \ roc / roc \ roc 1 |OQ 

WI V (L v2dX ) dX ~J VY (L y2dX ) dX ” /o y2U dX + 2 V 1o = ° (V - 33) 

The last term on the left is zero, and from continuity, 

roo / roo \ roo ^ IOO roo „ 

— J Vy IJ v 2 dx'\dx = U J v 2 dx\^-\- j v 2 Udx 


Thus, 



5 

ay 


f »(/>*') 


dx = 0, 


(V.34) 
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since U{ 0) = t>(0) = 0. 

Let 

F = J v (^J v 2 dx'^j dx (V.35) 

Then = 0, and F = constant, independent of the distance Y from the free surface. 
However, realistically, there may be some interference from the velocity leaving the 
surface and coming from the corner. Glauert states that physically, the “flux of the 
exterior momentum flux is constant, but this is hardly a familiar concept.” [Ref. 36]. 

It is worth pausing here to find a rough estimate of the magnitude of the flux of 
the exterior momentum flux, F. With a knowledge of the conditions in the impinging 
free jet, based on the entry conditions from the corner, start with 


F = Io° V 0f° ^ dx ') dx ~ V * fo° V (/°° w dx> ) dx ’ (V.36) 

where v* is a typical jet velocity. Write 


UM-HM 

(V.37) 

So, 


1 r°° d / r°° \ 2 

-aid vix ) dx 

(V.38) 

However, 




which concludes that 


F *\v-{j~ vdx f 

(V.39) 

This means that 


F = ^(typical velocity) x (mass flux) 2 
z 

(V.40) 


in the plane jet case. Neither the mass flux nor the magnitude of the fluid velocity 
will vary much when a plane jet is deflected off of the vertical wall, so (V.40) is a 
good estimate of the constant flux of exterior momentum flux, F. A more accurate 
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calculation of F will be given later. Keep in mind that this approximation for F has 
a singularity at the origin, where Y = 0. 

The similarity form of the velocity distribution depends only on F and is given 
Glauert’s paper [Ref. 36] as 

/"' + //" + 2/' 2 = 0 (V.41) 


with boundary conditions /(0) = /'(0) = 0 and /(oo) = 0. (As an aside, Batchelor 
refers to a steady narrow two-dimensional jet of fluid adjoining a plane rigid wall and 
offers 4/"' + //" + 2 f 2 = 0 as the differential equation [Ref. 37].) 

Start with Glauert’s equation, (V.41). Multiply the equation by /. Integrate 
the first term by parts: 

/ //'" = //" - Jf'f" = //" - \s a 

The second and third terms in the velocity distribution form a perfect differential, 

/ (ff" + vs' 2 ) = s 2 r 


Thus, 

//" ~\f + ff = const = 0, (V.42) 

since f'{oo) = 0. 

Multiply now by /“V 2 : 

r i/2 r - \r z, 2 f + f / 2 f =o. (v.43) 

Integrate again to obtain (this time, the first two terms form a perfect differential) 

f~ 1/2 f + |/ 3/2 = constant (V.44) 

o 

If fo(v) is a solution of the original third-order differential equation, then so 
is fi(rj) = Afo(Arj) for any constant A, and /i satisfies the boundary conditions if /o 
does. According to Glauert [Ref. 36], without loss of generality, select the solution 
with /(oo) = 1 and take the above constant to have the value of |. 
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Now, write / = g 2 . Then f = 2 gg' and the differential equation in / becomes 

= (V.45) 

which gives as a solution 

"= ln (^r)+^ ct “(^)- (V-46) 



Figure 7. Velocity Profile from Glauert’s Method 

A plot of the variation of Glauert’s velocity profile (/'(tj) vs 77 ) is given in 
Figure 7. This profile is similar to the profile determined using Timman’s approach 
and to the profile obtained using the numerical data from Canright [Ref. lj. 

The velocity and position equations may be 'written as 


(§f/» 

(V.47) 

( 5F y /4 * 

V32W 

(V.48) 
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and 



Recalling that rj = the above equation can be used to give an expression for the 
boundary-layer thickness, A: 


A = 


5 F \" 1/4 
32 Y 3 J 


(V.49) 


There is a drawback to Glauert’s method. The momentum flux has a singular¬ 
ity when Y = 0, which will not allow the velocity into the corner (along the surface 
at X = 0) to match the velocity out of the comer (along the wall at Y = 0). This 
means that the surface velocity is based on one ratio of the temperature gradient 
(u oc Tx 3 ), while Glauert’s suggested matching (V.40) bases the wall velocity on 
the inverse ratio (u a 7^ 1/3 ). Thus, along the surface, as the temperature profile is 
compressed, the velocity into the corner increases. However, along the wall, as the 
isotherms are compressed, the velocity decreases, which is not reasonable. 

As an alternative matching condition, to eliminate the singularity of the simi¬ 
larity solution at the origin, calculate F when Y is a distance equal to the boundary- 
layer thickness, A. This will give a strength value for F which is relatively equal to 
the value at the comer, but the singularity is now avoided. Let 


Solving for Y yields 


5 F 


(V.51) 


Now, matching the magnitude of the surface velocity expression into the corner to 
the maxi mum wall velocity expression out of the corner, we find 


f (ff= 2 - 5/3 ( 


/5 F\ 1/2 
\2Y ) 

(V.52) 

1/3 

1 

(V.53) 


which reduces to 


This mat ching condition for F embodies the positive feedback mechanism of the cold 
corner. 
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C. CORE FLOW REGION 

I NEVER COME ACROSS ONE OF LAPLACE’S “THUS IT PLAINLY APPEARS' WITHOUT FEEL¬ 
ING SURE THAT I HAVE HOURS OF HARD WORK BEFORE ME TO FILL UP THE CHASM AND 
FIND OUT AND SHOW HOW IT PLAINLY APPEARS. — N. Bowditch 


1. The Green’s Function 

Recall that the core velocity U is treated as incompressible irrotational flow. 
As explained earlier, a stream function ip is introduced such that V 2 ip = 0. Define 
a Green’s function G(x,x o) for Laplace’s equation. Thus, V 2 G(x, xo) = S(x — xq ), 
where x 0 is an arbitrary point in the domain. The solution can then be represented 
in integral form. First, however, the boundary conditions in the core flow region 
must be prescribed in terms of ip. The core flow must match the entrainment into 
the boundary layers, that is, the normal velocities. Along the wall boundary layer, 
U(x —* oo, Y) = ipy{X —» 0, Y). Along the surface boundary layer, V(X,y —> oo) = 
—ipx(X, Y —> 0). These conditions can be simplified to Dirichlet conditions, matching 
the stream function ip at the outer edges of the boundary layers. 

The Green’s function is determined using the method of images. Start with 
the solution in a quarter plane and use symmetry in the other quadrants to match 
the boundary conditions. There are four quadrants with logarithmic singularities. 
The defining problem for the Green’s function, V 2 G(x, x 0 ) = 8(x — x 0 ), satisfies the 
corresponding homogeneous boundary conditions, G(x, 0; xo, yo) = G( 0, y; x 0 , y 0 ) = 0. 
Here, the notation that G(x,x o) = G(x,y,x 0 ,yo) is used. The quarter-infinite space 
(x > 0, y > 0) has no sources except at a concentrated source at x = xo. 

Using the method of images in the core flow region, the Green’s function for 
Laplace’s Equation on a quarter-plane (X > 0, Y > 0) is 

G(x,x 0 ) = i-{in[(x-x 0 ) 2 + (y-y 0 ) 2 ]-in[(jf-xo ) 2 + (y+n) 2 ] 

-In [(X + X 0 ) 2 + (Y- V 0 ) 2 ] + In [(X + X 0 f + (Y + V 0 ) 2 ]} 
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- 1 ■ f t(£ - *o) 2 + <x - r ») 2 ) k*+ x 0 f+< y+ r 0 ) 2 n 

■fcr U(x-Xo) 2 + (v' + yo) 2 ][(x+Xo) 2 + (r-r 0 ) 2 ]/’ 1 ' ’ 

To check the boundary conditions, start with the Green’s function and set X 
or Y, respectively, equal to zero. Thus it plainly appears that G = C? = 0. 
Since the boundary condition for the Green’s function is the Dirichlet kind at Y — 0, 
and at X = 0, then a positive image source must be used for Y <0 and X < 0, and 
negative sources must be used for X > 0, Y <0 and for X < 0, Y > 0. In this way, 
equal but opposite sources are located at the symmetry locations of Xo in the four 
quadrants of the plane, and there is no flow along Y = 0 and X = 0, as desired. 

To solve Laplace’s equation, multiply the Green’s function identity by ip, giving 

[v 2 G(X,Xo) = 6(X-Xo)]. 

Multiply Laplace’s equation by G, yielding 

G [VV = 0]. 

Subtracting and integrating over the domain A gives 

f f V 2 G - G VV) dA = (Xo) (V.55) 

Using Green’s second identity, Laplace’s equation in two dimensions becomes: 

I J a V 2 G - G Vfy) dA = jf (GVxp -ijfVG) •nds = 'ip (X 0 ) , (V.56) 

where h is the normal to the boundary (inward normal) and s is the arc length along 
the boundary. Consider a large quarter-circle (in the limit as the radius tends to 
infinity). The inward normal for the top (boundary layer along the free surface) is j 
and the inward normal for the side (boundary layer along the wall) is i. It will be 
shown that the contribution at oo tends to vanish if ?/> —> 0 sufficiently fast. 
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(V.57) 


+ lim <f(GXiP-ipVG)-nds 

|X|—oo Js 

The last term on the right is the contribution as the radius of the quarter circle (or 
|X|) tends to infinity. It will vanish if ip —► 0 sufficiently fast, or if ^ ~ To show 
this, integrate in polar coordinates, for which ds = rdd (r = |X|): 


hm f (GV^ -ipXG) rdO = 0. 

Substituting G and evaluating § d6 = 2ir yields the necessary condition 

( dip \ 
hm r lnr —- ip) =0. 


(V.58) 


dr 


(V.59) 


If ^ 0 sufficiently fast, or if ^ ~ then the contributions of this “boundary” 

\X\ 

term wifi vanish. 

Hence, 

dip 




dG 

dX 


dY + 


x=o 


r{ a w-*w) L" < v - 6o) 


Since G 


x=o 


= G\= 0 , 


y=o 


* w -r h n) l <nr+ r{-*w) 


\x=o 

and V = -#& 

G = 0 at X = Y = 0, and 


dX 


Y =0 


(V.61) 


Recall that V = U 

dY x=o 


dX 


y=o 


. Since the boundary conditions give 


= 0, we only need to find dG X x °) which is 
x=o ’ J 9Y y=o 


(V.62) 


dG(X,X 0 ) 

_ 1 

F 0 Y 0 

dY 

y=o — t r 

[(x + x 0 ) 2 + r 0 2 (x - x Q f + r 0 2 J 


Thus, for the core flow due to the surface boundary layer, insert the mass flux into 
the corner for yielding 

* = -/%W*,o) 3G( * ,J?o) i 

Jo 


dY 


y=o 


dX, 


(V.63) 


As stated previously, one advantage of the Green’s function is that the flow 
can still be represented over the entire quarter plane, so that there is no artificial 
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recirculation due to imposed artificial boundaries. Thus, only the boundary condi¬ 
tions at the boundary layer/core region interfaces are required. The Green’s function 
approach does not require flow boundary conditions at the computational domain’s 
boundaries. 


2. Mass Flux Calculations 

In order to obtain a solution to Laplace’s Equation in the core flow region, the 
boundary conditions from the surface and wall boundary layers must be addressed. 
Recall that along the surface boundary layer, the velocity is modeled as 

u = ^ e -1/3 (-T) 1/3 (Tx) 1 ' 3 erfcfr). (V.64) 


As explained previously, the horizontal velocity component u decays to zero as X and 
rj increase. 

To determine the mass flux along the surface layer into the cold corner, sum the 
velocity component as (the correct scalings for the velocity and the boundary-layer 
thickness are inserted as appropriate) 


^surface 



(—T) 1 / 3 (T x ) 1 / 3 erfc( ?7 ) 


dy 


= ~ Jo 1/3 ( _T ’) 1/3 ( T ^) 1/3erfc ( 7 ?)] 6dr l 




= -^- 2/3 (-r) 2/3 enr 1/3 f 


(V.65) 


e~ 213 (—T) 2/3 (T x )~ 1/3 

A 
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When X = 0, the stream function (mass flux) into the corner can be calculated as 
\ s -2 / 3 (—To) 2 / 3 (T Xo )~ 1/3 . This value will be compared to the mass flux out of the 
corner shortly. 

To determine the mass flux down the wall, sum the vertical velocity compo¬ 
nent, v. Recall that 



where 

x / 5 F \ 1/4 
^ = A = v32yv *• 

Thus, the stream function (mass flux with appropriate scalings for the velocity and 
boundary-layer thickness) is 



= -(5F) > ' 4 (2) 3 ' ,4 r 1 / 4 /{>?)|” (V.66) 

We have shown that /(0) = 0 and lim^oo f(r}) = 1, so that the mass flux down the 
wall is 

^waii = -(40 F7) 1/4 (V.67) 


D. POTENTIAL FLOW DUE TO A WALL JET 

The next test for our problem is to check the potential flow due to a wall 
jet in the quarter-plane. A jet is a flow in which the width, or cross-stream scale, is 
considerably smaller than the downstream scale. For a wall jet, this type of flow occurs 
along a solid boundary. In Glauert’s discussion of the wall jet, he postulated that the 
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entire flow of the wall jet cannot conform to one similarity solution over the entire 
domain [Ref. 36]. Glauert divided the flow region into inner and outer portions, 
treating the two portions separately. The dividing line for the two portions was 
the point of maximum velocity. Additionally, as we have demonstrated previously, 
Glauert also showed that the “flux of exterior momentum flux” across the viscous 
boundary layer is constant [Ref. 36]. The plane wall jet is created when a jet flow 
is discharged in a tangential direction along the plane surface; in this case, the wall. 
Plotkin [Ref. 40] studied a higher-order correction to Glauert’s solution, providing a 
second-order solution including displacement. The first-order inner viscous solution 
induces a correction in the outer inviscid core flow solution. Plotkin matched the inner 
and outer solutions and developed a solution which satifies the boundary conditions 
in a semi- infini te plane, where ip(x, 0) = 4x 1//4 for x > 0, and ^ = 0 for x < 0. His 
solution satisfies Laplace’s equation and is 


ip(x, y) =4 ^R.(x + iy) 1 ^ — S(z + iy) 1 / 4 ] (V.68) 


where K and S denote the real and imaginary parts, respectively, for 0 < arg(x-Hy) < 
27 r. 


In the present problem, Plotldn’s solution must be modified to match the 
boundary conditions in the quarter-plane, namely, that ^(0, Y) — —{AOFY) 1 ^ 4 and 
when Y = 0 and X > 0, V’ = 0. The angle is |, and after taking the fourth-root, the 
argument for the complex variable becomes |. Taking into account that the jet flows 
down the wall (in the y-direction), the potential solution becomes 


1>(X,Y) = —(40 F) 1/4 


K{Y + iX) 1/A - 


V2 + vg a(y +, X ), /4 

y 2 — y/2 


(V.69) 


which simplifies to 


1>(X, Y ) = -(40F) 1/4 [K{Y + iX) 1/4 ~(l + V2) 3(Y + iX) l/A ] . (V.70) 


This also satisfies Laplace’s equation in the quarter-plane. This expression for ip(X, Y ) 
is used to determine the core velocity field due to the wall boundary layer. 
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VI. NUMERICAL RESULTS 


REFINEMENT OF THE METHOD ALSO REQUIRES A MORE REFINED USER. — Federico Paris 
and Jose Canas 


The PURPOSE OF COMPUTATION IS INSIGHT, NOT NUMBERS. — Richard W. Hamming 


A. SOLVING LAPLACE’S EQUATION 

Before solving Laplace’s equation in the core-flow region, the numerical code to 
solve Laplace’s Equation was tested on a simple problem with a known solution, the 
point vortex. The two-dimensional solenoidal flow in the core flow region associated 
with a straight line vortex may be described in terms of a stream function as 


ip = —^loga, (VL1) 

where k is the strength of circulation and a = \J{x — £ 0 ) 2 + (y — Vo) 2 - In a flow field 
which is entirely two-dimensional, the appropriate term for the singularity is point 
vortex. The line vortex, a line in space that has a given circulation, or vorticity, 
becomes a point vortex when a single plane cuts through it. For simplicity, set k = 1 
and find 0, the harmonic conjugate of ip. Since the two scalar functions <p(x, y) and 
ip(x, y) both satisfy Laplace’s equation, they also provide alternative specifications of 
a core flow vector U which is both irrotational and solenoidal [Ref. 37]. Thus, given 
ip above, and the fact that 


dip 

dy 


d(p 

dx 


(p is found to be 


K, 

2tx 


,-i 


&P =y _ d<P 

dx dy ’ 

(VL2) 

6 ) 

(VI.3) 
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To find the velocity anywhere in the fluid, find the distance and direction to 
the point vortex. The direction of the fluid is always a circle surrounding the point. 
The magnitude of the velocity is given by 

M - 2ST (VL4) 

The horizontal and vertical components of velocity, evaluated at x = 0 and y = 0, 
respectively, axe 

u = — 1 v = .1 

dx x=o y ’ dy j/=° x 

The Laplace’s equation solver was rim on this point vortex problem, treating 
the point vortex at the origin and evaluating the <j> distribution in a quarter-plane with 
a uniform grid. A trapezoidal scheme for integration is used. The output matched the 
behavior of the known solution, as seen in Figure 8. The plot of the tan _1 (y/x) gives 
straight fines through the origin. The numerical code integrates to regions beyond 
the area of interest (due to the Green’s function), and then the area of interest - the 
comer region - can be analyzed in detail. This was useful to test the code, although 
the stream function ip will be used instead of the velocity potential <p. 



B. 


SOLVING THE HEAT EQUATION 

The full two-dimensional time-dependent problem can be expanded as 


far dT v dT\ d 2 T d 2 T 

l. at + u ax + v ay) ~ ax 2 + ar 2 ' 


(VI. 6) 


The boundary conditions for the quarter-plane are that when y = 0, T y = 0, and 
when x = 0, T = — 1. Artificial boundaries were introduced to modify the quarter- 
plane problem into a rectangle problem. Since we are concerned with the activity in 
the corner, the domain was chosen large enough that the temperature gradients and 
velocity values decayed to zero far from the corner (as X, Y —> 0, VT, U, V —>■ 0). 
Specifically, on the right boundary, T — 0 for the hot incoming fluid, and on the 
bottom boundary, Ty = 0 to minimi ze downstream influence. Along the surface layer, 
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Figure 8. Contour Plot of Numerically-Solved Phi for Point Vortex 


the heat flux Ty condition (IV. 17) is incorporated into the scheme. It is assumed that 
the time-dependent solution will converge to a steady-state solution. Hence, the initial 
conditions can be chosen arbitrarily. 

Equation (VI.6) was solved using an alternating direction implicit (ADI) method 
as developed by Polezhaev in 1967 [Ref. 39]. This method treats each time step as 
two half steps. The above equation can be substituted into the scheme as: 


tj. . _ p 

Ul ' 3 2 AX dx , 


STEP 1: 

, A t ( 

1+ t( 

STEP 2: 

1 At (ir $Y c 2 

1 + T \ Vi ’’ 2aF “ Sy 


rp* _ 

““ 


1 - 


At (y. 6y p" 

V l ’ 3 2A Y dy , 


rpn+1 


'-IK 


rpn 

1 h3 


2 AX ~ 6x , 


h3 


(VI-7) 

(VI-8) 
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where the central-difference operators <5 and <5 2 are defined as: 

8x Ti,j = Ti + ij — Ti-ij 

8y T it j = T it j +1 — Tij-i 


(VI-9) 


= 


rnn o T^n i ^Pn 


(VI. 10) 


x iJ (AX) 2 

rpn _ o rnn i r PJi 

= wi (A y ^ (vi-io) 

As a result of splitting the time step in this algorithm, only tridiagonal systems 
of linear algebraic equations must be solved (at each half-step). This method is 
first-order accurate in time but second-order in space, with a truncation error of 
0[At, (Ax ) 2 , (Ay) 2 } and is unconditionally stable for the linear (no convection) case. 

The scheme was coded using the C programming language with a uniform 
(Ax = Ay) rectangular grid. The problem’s boundary conditions as stated earlier 
were incorporated into the scheme. A copy of this code is attached in Appendix E. 
To set the boundary conditions at each of the four sides of the rectangle, some of 
the coefficients for and T™* 1 had to be manipulated. These manipulations 

ensured adherence to the boundary conditions without altering the solution method 
inside the boundaries. For example, second-order difference equations for the first 
and second spatial derivatives were modified at the surface and at the bottom of the 


A subroutine tridiag. c is called in the execution, which solves the tridiagonal 
set of equations Ax = b for the unknown left-hand-side of each equation in STEP 1 
and STEP 2 of the ADI method. In this case, the matrix A is the coefficient matrix 
involving velocity and time-derivative terms, b is the vector consisting of the right- 
hand-side of each equation in STEP 1 and STEP 2 (the known temperature), and x 
is the vector of the unknown temperature. 

The method was tested against known problems involving mixed boundary 
conditions and various scenarios of velocity, with and without convective terms, and 
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after iteration, a steady-state was reached, i.e., the solution converged. The ADI 
method is unconditionally stable but may converge slowly. An alternate technique is 
to update each row or column in the step temperature as it is solved. This leads to a 
faster convergence; however, it is not clear if this modified technique is unconditionally 
stable. 


The first sample problems that the ADI method were tested against involved 
zero velocity (|| = V 2 T - pure diffusion, no convection). The solutions were linear in 
x : T(x, y) = where L is the distance to the boundary in the rr-direction. Other 
sample problems involved a non-zero velocity in one direction + u • VT = V 2 T 
- diffusion and convection). In the latter case, one example used involved a uni t 
square domain, with Dirichlet conditions (T = 0) on the left and right boundaries 
and Neumann conditions on the bottom ( T y = 0) and surface: 


f(x) = e 1 / 2 sin(7rx) sinh 


/ \/l + 47r 2 ' 


V 


The velocity field was (U, V ) = (0, —1) in this example. The solution to the problem 
is 

(Vl+4*-2-l)(y-l) (V'l+4Tr2+l)(y— 1) * 

e 2 e 2 

Vl + 47T 2 — 1 \/l +47T 2 + 1 



The ADI algorithm converged to four-decimal accuracy in 49 iterations, yielding the 
solution shown in Figure 9. The grid had a uniform spacing of 0.05. The time step 
per iteration was 0.002. 


C. THE ADI METHOD 

The ADI method was used to solve the two-dimensional unsteady heat equa¬ 
tion, beginning with an initial temperature field which was exponentially decaying 
on the surface and zero elsewhere. This temperature field was then iterated in two 
loops: the inner loop had a fixed number of iterations which was simply the AP T 
algorithm marching in time; the outer loop updated the velocity profile with another 
fixed number of iterations. Since the velocity field is based upon the temperature 
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After 50 Iterations - ADI method 



Figure 9. Contour Plot of Numerically-Solved Temperature for Test Case 

and its derivatives, the interim temperature solution (after each set of inner loop it¬ 
erations) was used to calculate the velocity, which was then fed into the ADI solver 
for another series of inner loop iterations, yielding another, more accurate, tempera¬ 
ture solution. This process was repeated a number of times equal to number of the 
outer loop iterations. Once the difference between successive temperature solutions 
fell below a certain tolerance (in this case, 0.00001), steady-state was reached. Figure 
10 shows a typical steady-state temperature solution. The spatial step size is 0.01 in 
both directions. 

The thermal field is compressed against the wall by the velocity into the corner. 
The problem used 20 iterations in each inner loop, and it converged in only 12 outer 
loop iterations. A closer view of the corner is shown in Figure 11, which also compares 
our approach (on the left) with Canright’s numerical data (on the right). The com¬ 
parison is encouraging, particularly since the model has no adjustable parameters. 
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Steady-State Temperature Solution 



Figure 10. Steady State Temperature Contour Plot 


To get a better indication of our method, surface temperatures are plotted, in 
Figure 12. There is excellent agreement of the temperature gradient from the corner 
into the core flow region, well beyond the boundary-layer thickness. 


D. UNIFORM VELOCITY 

Following the boundary-layer matching theory outlined in Bender & Orszag, 
[Ref. 41], the uniform velocity for the entire domain can be written as 


^uniform — dinner d~ 


outer 


^match 


(Vl.ll) 


The inner solution is composed of the velocity components in the two boundary layers, 

tanner = ^surface d" ^Avall, (VI. 12) 


59 




Figure 11. Steady State Temperature Contour Plot - Comparison with Canright’s 
Data 


while the outer solution is made up of the two components of the core velocity, due 
to the similarity solution and the Green’s function: 


U c 


outer 


— Ugreen H” C^sim U a 


(VI. 13) 


In the above expressions, £/ sur face refers to the surface boundary layer velocity 
component obtained using the modified Timman approach, due to the thermal gradi¬ 
ent along the surface, T x . t/ wa n refers to the wall boundary layer velocity component 
obtained using Glauert’s method, fed from the surface mass flux. The core velocity 
component Ugreen is a velocity field obtained using the Green’s function (from solving 
Laplace’s equation in the core region) and due to the entrainment into the surface 
layer, and the core component U s i m is the velocity field due to the similarity solution 
of the potential from the wall jet (due to the Plotkin modification). Determine the 
uniform flow in terms of the stream function, ■0. 

Consider first the surface velocity vector. As shown earlier, 

tWfeo = -^£-‘' 3 (-r) 1/3 (Tx) ,/3 erfc(i;). (VI. 14) 
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Steady-State Surface Temperature 



Figure 12. Steady State Surface Temperature Plot 


Let g(rj) = e 1)2 ds = & erfc(7?), where v = y/ -pjv . Define G (r?) = g(s) ds = 


1/3 


f (1 — e n ) + yg(v)- Now, the stream function may be written as 


ry (-T) 2 r Y 

Surface = / lldj/ = - G{rj) = / U dY’ 

JO £' lx JO 


21 V3 


(VI.15) 


The scalings are V = — V’X) ^ = u/P, and Y = Py. 

Next, consider the wall velocity component. The stream function may be 
written as 

^wai. = - j*v dx' = - J* V dX' = -(^py) 1 / 4 /(r?), (VI.16) 


where 


v = 


T] = 


r5FlV2 rw 


2yj 

5 F 1 1/4 


32y 3 


X 
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The scalings are U — ipy, V = v/P, and X = Px. Here F is found from (V.53) to 
match the surface flow to the wall jet. 

For the similarity component to the core velocity, ?7 s i m , which is due to the 
wall jet, the stream function may be written as 


^sim = - [40F] 1/4 [U(Y + iX) 1/A - (l + y/2) 9f(y + iX) 1 ' 4 } , (VI. 17) 

with the two boundary conditions that ^(0, Y) = — [40F] 1 / 4 Y 1 ^ 4 and ip(X, 0) = 0. 

For the potential component to the core velocity using the Green’s function, 
U green, the stream function may be written as 


i’l 


green 


. r°° 1 

= -Jo 2 


f(-r) 2 l 

1/3 dG(X,X o) 

U 2 Tx\ 

dy 


y—0 


dX , 


(VI. 18) 


where the Green’s function is 


G(X,X o) = ^{to [(* - *o) 2 + (y- Vo ) 2 ] - to [(x - x 0 ) 2 + (y + yo ) 2 ] 
-In [(x + x 0 ) 2 + (y- yo) 2 ] + to [(x + x 0 ) 2 + (y + yo) 2 ] | 




- so) 2 + (y - yof] [(g + ^o) 2 + (y + y 0 ) 2 ] 
o) 2 + (y + yo) 2 ] [(* + xo ) 2 + (y- yo) 2 ] 


and 


dGjXj C p)| 
dy 


• _ 1 

yo 

yo 

*y=o 7r 

(x + x 0 ) 2 + y$ 

(x-x 0 ) 2 + ^. 


(VI. 19) 


(VI.20) 


Finally, s umm in g all of the components yields the uniform velocity, which can 
be written as 


^uniform (X, Y) = Y) + 1pgreen(X, Y) + ^surface (V, Y/P) + 1p wa ,\\(X/P, Y) 

'0surface(-^> Oo) 1pwall(.&- > iY) (VI.21) 

The Prandtl number P is simply a parameter in the model that represents how 
thin the boundary layers are. The boundary-layer model should be independent of P. 
Using the uniform stream function, a streamline field can be plotted. This is shown 
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in Figure 13. The horizontal velocity component increases as the distance to the 
wall decreases, and the vertical wall jet can clearly be seen. The flow into the corner 
along the surface turns down the wall, due to the jet. The upper left figure shows the 
complete contour plot of the stream function in the domain. The upper right figure 
shows the core stream function, Y). The middle set of figures represent that 

portion of the stream function which is due to the wall components, ^ WSL n(X, Y) + 
Y)—i> waii(oo, Y) on the left and the similarity solution ip S im(X, Y) on the right. 
The bottom row of figures shows those portions of the uniform stream function which 
are due to the surface components Y) + ^Wfecepf, Y/P) - i> sw{ace (X, oo)) 

on the left and that due to the Green’s function (?/' green (X, T) on the right. A closer 
detail of the corner stream function dynamics is shown in Figure 14. 

As can be seen, the velocity field appears to be strong enough to compress 
the thermal field into the cold corner. The wall jet similarity solution dominates 
the Green’s function solution in the core flow region and overall in the steady-state 
solution. The few streamlines near the corner in the Green’s function contour plot 
suggest that h(x) = ^ reaches a local maximum along the surface, not far 

from the corner, and the impact of h(x) is strong enough to alter the velocity field in 
this view. However, overall, h(x) does not appear to change the expected behavior of 
the velocity field. 


E. UNIFORM TEMPERATURE 

In a similar fashion to finding the uniform velocity approximation, the unif orm 
temperature for the entire domain can be written as 


* uniform 


(VL22) 


where Tj nner is the temperature solution pertaining to the surface boundary layer, 
T out er is the numerically-obtained temperature solution from the ADI method (solving 
the unsteady two-dimensional heat equation with periodic velocity updates), and 
^match is the matching condition which satisfies the flux condition at the surface. 
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From Surface 


Steady-Stale Uniform Stream Function 


Steady-State Uniform Stream Function Due to Core Components 





Figure 13. Uniform Stream Function and its Components 
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Figure 14. Uniform Stream Function and its Components — Closer Detail 
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Recall that the dominant balance argument gives uTx — t yy as the surface heat 
equation (see Chapter IV). 

Integrating the surface heat equation through the boundary layer gives 

ty(X,y)-t y (X, 0) = FuT x dy' *T x {X, 0) f udy' 

J 0 i/0 

= T x (X,0)4>(X,y) = t,(.X,y) + T Y (X,0) (VI.23) 

As y —» oc, t y —• 0, and the heat flux condition in the core region becomes 

T y (X,0) = T x (X,0)4’(X,oo) (VI.24) 


Solving for t y yields 

roc 

t y (X,y)=Tx(X,0){1>(X,y)-MX,oc)] = -T x (X,0) udy' (VI.25) 

Jy 

Integrating this expression from y to oo gives 

t(X, oo) - t(X, y) = -T x (X, 0) J™ (jT u dy ") dy' (VI.26) 


or 


t(X,y) = T x (X, 0) r r u d/ dy' (VI.27) 

Jy Jy f 

Prom Timman’s method, the horizontal component of velocity is u = — j X- T j r -v j ^ 
where rj = y/6, 6 = , G(y) = $g(0 d£, and g(y) = J*~ e~ s2 ds = ^eric(y). 


Thus, 


t y — Tx V’oo) 
(- T)T x l2/3 


V5F 


roc 

/ erfc(r/) dy' 

Jr) 


(-T)Tx 


I 2/3 


ie 7,2 -y g(y) 


(VI.28) 


g(v )> 
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and the temperature t may be written as 


t = 


t x r (<i> - m <v 

= (17) [( 1+2 ^)s( 7 ?)- ? ? e " V ] 


(VI.29) 


Uniform T*mp*r«tur» Solution 



As can be seen in Figure 15, the uniform temperature solution (dotted contour 
lines) matches the steady-state solution (solid contour lines) far from the surface (as 
77 gets large), but there is a slight difference when 77 is small (close to the surface), 
due to the G(0) term. The uniform temperature prediction bends toward the surface 
to intersect it perpendicularly, satisfying the T y = 0 boundary condition. 


F. VORTICITY 

The vorticity plays a major role in the cold corner, turning the flow from the 
surface down the wall. Our model has no vorticity in the core flow region; however, 
the numerical data shows some vorticity there, so the flow fields look different. T his is 
due to the hypothesis that the vorticity that is created far upstream from the corner 
(along the surface) diffuses into the slow flow and gets transported to the region near 
the cold corner and not into the boundary layer regions. Figure 16 shows a contour 
plot of the stream function (dotted contour lines) and vorticity (solid contour lines) in 
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the cold corner, due to Canright’s numerical data. Most of the vorticity is confined to 
thin regions along the surface and wall (the boundary layers), but some vorticity is in 
the core flow region (where the streamlines are sufficiently far apart to designate slow 
flow). The main result is that core vorticity is not essential to the feedback process, 
based on the comparison of the thermal fields of the model and numerical results. 

Vorticity and Stream Function Contour Plot 



Figure 16. Vorticity and Stream Function Contour Plot (due to Canright) 

To quantify Figure 16, far from the corner (x ~ 3), the surface value of vorticity 
is O(10 3 ) smaller than than surface vorticity value at the corner (x —> 0). Recall that 
a; = 0atr = y = 0. Even at a distance of <5 from the corner (from both the surface 
and the wall), the vorticity value is O(10 2 ) smaller than the peak value in the corner. 
Thus, the vorticity outside the boundary layers is non-zero but insignificant compared 
to the boundary-layer values and is not essential to our model. 
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VII. 


CONCLUSIONS AND DISCUSSION 


Complex Happens. — David Canright 


The thermocapillary stress condition at the free surface gives a starting point 
for the scalings and dominant balance. We have obtained relations for the thermal 
length scale l, the viscous thickness <5, the surface velocity u, and the core velocity U 
scales in the cold corner, in terms of the dimensionless parameters M and P. These 
scalings are 


6 ~ M~ l P~\ 
l ~ M-'P- 2 , 
u ~ P, 

U ~ P 2 . 


These scalings are consistent with the numerical results obtained by Canright in his 
1994 work [Ref. 1 ]. 

The temperature gradient, Tx, on the surface drives the strong surface flow, 
similar to a surface jet. In the viscous surface boundary layer, a modification to 
Timman’s method was employed which allowed the surface velocity to be modeled as 

u = e ~ 1/3 (—X 1 ) 1 / 3 (T x ) 1 / 3 erfc( 77 ), (VII. 1) 


where e is a constant and 77 — y/8. The boundary-layer thickness 6 in terms of the 
surface temperature is 



(VII. 2 ) 


The surface flow then feeds into the corner and into the wall flow, similar to a 
wall jet. In the viscous wall boundary layer, a modification to Glauert’s method was 
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employed which allowed the wall velocity to be modeled as 

/5 F\ 1/2 

v ={w) m 


(VII.3) 


where 


' 5 f \ 1/4 

,32 Y 3 ) '' 


(VII.4) 


and F, the flux of exterior momentum flux, is found from matching the surface velocity 


16 y/n (T x 


m ■ 


(VII.5) 


Using these expressions for the velocity components, the stream functions 
were calculated. The strong flow in the viscous boundary layers induces a relatively 
weak irrotational flow in the core region. Two approaches were employed to find the 
core velocities: a Laplace’s equation solver using a Green’s function, and a complex- 
variables solution modeled after a self-similar plane wall jet. These velocities were 
then substituted into the unsteady heat equation and an alternating-direction implicit 
(ADI) numerical scheme was employed to solve for the temperature. The weak core 
flow then contains a thermal gradient which was used to determine the Tx on the 
surface. Once found, the surface temperature gradient was used to update the surface 
velocity and the process was iterated until steady state was reached. 

When the flow along the wall is modeled as a wall jet, a modification to 
Plotkin’s second-order approach was used to match our boundary conditions in the 
quarter-plane. Thus, the potential solution becomes 


ip(X, Y) = -4 [lZ(Y + iX) 1/4 -(l + V2) %(Y + iX ) 1 ' 4 ]. (VII.6) 

The velocity components may be found by U = and V = — J^. 

It is fo un d that the flow into the cold comer is strong enough to contain 
the thermal field. The isotherms are compressed along the wall after steady-state is 
reached. The temperature profile obtained from the ADI method compares favorably 
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with Canright’s numerical data from 1994. Considering the simplicity of our model, 
the results are very encouraging. The thermal field results are self-consistent. 

The prediction model of the uniform velocity vector field is encouraging, in 
that it matches the inner, boundary-layer flow to the outer, core flow, and subtracts 
the matching flow elements. The basic model applies for the whole range of Marangoni 
and Prandtl numbers in the convective-inertial regime. The unif orm flow picture will 
change with a changing Prandtl number, reflecting the relative thickness of the viscous 
boundary layers, but the underlying model does not change. In addition, because this 
model seems to capture the cold corner feedback process, these results could be useful 
as a local solution in a larger numerical problem, without the need to numerically 
resolve the viscous boundary layers. 
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VIII. FUTURE WORK 


We ask, with the extension of mathematical knowledge will it not finally 

BECOME IMPOSSIBLE FOR THE SINGLE INVESTIGATOR TO EMBRACE ALL THE DEPARTMENTS 
OF THIS KNOWLEDGE? — David Hilbert 


There are a few limitations to this work. The grid used in the numerical 
approximations is uniform. To get a more accurate description of the velocity and 
temperature fields, a non-uniform grid should be employed, which is finer near the 
corner and more coarse away from the corner, as in Canright’s 1994 work. Even so, 
the current data shows good correlation with the previous work. A nonuniform grid 
would require a more complicated algorithm but may run faster for a given resolution 
in the corner, especially since running of the current ADI codes with updates required 
hours for several thousand iterations. 

More effort should be placed on the effects of vorticity in the corner. In the 
numerical data, the vorticity does not appear to be entirely confined to the boundary 
layers, yet the vorticity in the core flow region seems to have no significant impact 
on the temperature field solution. It appears that the core flow does not turn down 
the wall if vorticity is not present. The vorticity in the surface and wall boundary 
layers must be strong enough to force the nonlinear flow at the free surface down the 
rigid wall. At this point, we can estimate the magnitude of the vorticity, but more 
research is needed to determine how to combine the vorticity into the potential flow. A 
meaningful approach to modeling the vorticity throughout the domain would provide 
a good understanding of the rotation mechanism in the flow. T his was addressed in 
Chapter VI. 

This research simplified the physical problem. The geometry was revised, some 
material properties were ignored, and the free surface was assumed flat due to the 
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amount of surface tension. If a geometry change occurs, that should be investigated 
in future research. 

Also, higher-order modifications to Timman’s method should be attempted. 
This would provide more accuracy near the free surface. The current model predicts 
the decay of the velocity through the boundary layer fairly well, until the interface 
with the core region is approached, and then the Timman modification decays faster. 
If a four-parameter model could be developed, it might better depict the velocity 
profile using the thermocapillary condition at the surface. 

As more experiments are conducted in the low-gravity environment, such as 
those using the NASA space shuttles by Kamotani and Ostrach and others, a valida¬ 
tion of the current work using experimental techniques should be attempted. 
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APPENDIX A. TIMMAN’S METHOD ALONG 

THE SURFACE 

Even though the one-parameter model is very encouraging, the goal is to try 
to derive a more accurate velocity profile than the simplest approach in Chapter V, 
since the surface tension and thermocapillary condition at the surface is driving the 
flow. Now, suppose the velocity profile for u(y) (where the A dependence is implicit) 
is 

u(y) = f(tj) = jf° e - * 2 (a(A) + c(X)t 2 ) dt - b{X)e~ Tt \ (A.l) 

This introduces a damping term into the velocity profile. Is there a more accurate 
expression for the boundary-layer thickness? At the surface, y = 7 ? = 0, which implies 
that 

ti(0) = /(0) = e~ t2 (a(X) + c(X)t 2 ) dt - b(X) 

- ^«P0 + ^o(X) - b(X) (A.2) 

In the core flow region, as y —*■ 00 , 77 —> 00 , and /( 00 ) = 0 . 

Differentiating f'(rj ) = e^ 2 [-a(X) - c(X)r) 2 + 2b(X)r]}. Evaluating 

this at the surface gives /'( 0 ) = -a, which yields a(X) = -6(X)T x , just as in 
Chapter V (when 6 (A) = c(X) = 0). 

Differentiating further yields f"(r)), which when evaluated at the surface gives 
/"(0) = 26(A). Since f(rj) = u, f'(rj) = 6(X)u y , and f"{r}) = 6 2 (X)uyy. Thus, 
6 (A) = 1 6 2 (A) u( 0) ttx( 0 ) at the surface. 

Taking the third derivative of /(??), /"'(r?) is a very lengthy expression. u yyy 
can be rewritten as ~ (u yy ), which after substitution into the momentum equation 
becomes |(uu, 4- V u y ). This can be simplified using the continuity equation and 
the fact that r( 0 ) = 0 . Evaluation at the surface yields /"'( 0 ) = 6 Z u(0) T X x- Since 
f'"( 0) = 2a(X) — 2c(A), we can solve for c(A). 
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The proposed velocity profile at the surface is therefore 

«•(#) = f(v) = -j°° [s(X)T x + { ST X + ±6 s u(0)Txx}t 2 ] e- ,2 dt-±6 2 (X)u(0)u x (0)e-” 2 . 

(A-3) 

Hence, a first-order, nonlinear, nonhomogeneous ordinary differential equation for the 
velocity can be developed in terms of the boundary-layer thickness and the tempera¬ 
ture gradient as (let u(0) = uq): 

«o«o + «0 (| + = -- W Tx (A - 4) 

where u' 0 = ^u(O). This gives one equation with two unknowns, uq(X) and S(X). 

We also have <5 2 = —T = 6 / 0 °° u 2 d rj. This type of ordinary differential equation (A.4) 
can be known as an Abel Equation of the Second Kind [Ref. 33]. Finding a solution 
to this type of equation is available only for certain cases, such as if the coefficients of 
uqu' q and uo fall into special categories outlined in [Ref. 34, 35]. Unfortunately, this 
particular ODE does not fall into one of those special cases. 

A second approach to modeling the velocity with a Timman-type method of 
two parameters is algebraic. Write the velocity profile as 

u = f{rj) = f ae~ t2 dt — be* 7 ?. 

JT) 

(Remember that the variables a and b are functions of X; however, for simplicity, 
the a(X) and b(X) will be replaced by a and b ). Using three derivatives evaluated 
at y = 0 will enable us to obtain an expression for u in terms of the boundary-layer 
thickness and then in only temperature terms. Recall that 

m =-» 

Taking derivatives of / with respect to 77 and using the chain rule since y = Srj, we 
obtain 

f(y) = e -7?2 (—a + 2by) = 8u y 

/'(0) = -a = 8(u y ) y=0 = 8T x (A.5) 
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Thus, a = — 6Tx • Continuing, 


f"(v) = e- r) \2b + 2a V -4br ] 2 ) = 8 2 Uyy 
/"(0) = 2b 

f'(v) = e- T > 2 [2a-8br l -2r ] (2b + 2a V -4b n 2 )}=6 3 u yyy 
0) = 2a — —26 Tx 

Recall that u yyy = ^ (u yy ) = (uux + Vu y ) = uTxx- Thus 

/'"(0) = 6 3 u(0)T xx = 5 3 (£. - bj T xx = « 3 
which, when equating the /'"(0) terms, finally yields 

= + (A.8) 


(A.6) 

(A.7) 


so that 


6 = 


2T X _ V^crp 

6 2 T XX 2 STx 


(A.9) 


and the velocity profile can now be modeled as 

2 T x 


u 


= - r STx e~ t2 dt- 
Jt] 


S*T XX - T STx 


o-*r 




8 Tx erfc(7/) 


2 T x y/n 


8 2 T 


8T- 


x 


XX 


o-r 


Recalhng (V.9), the boundary-layer thickness 8 can be solved as 

«(*) = 75^3- 

Jo u dr ! 

This yields a quadratic equation in <5 3 , which is: 


6 6 \t\ 




+ s 3 


(A. 10) 


(A.11) 


At the corner, with T, Tx, and Txx values taken from the numerical data 
at x = 0, the boundary-layer thickness 8(X) is 0(M~ 1 P ~ 1 ), as expected. As X 


77 


increases from the wall, 6(X) decreases in value, also as expected. However, after 
a certain value of X, still within the comer region, the solution to (A.12) results in 
complex roots. Unfortunately, this means that the algebraic two-parameter approach 
had to be abandoned. 
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APPENDIX B. CALCULATION OF e 


In the surface boundary-layer thickness equation, (V.22), the constant 

7T 

£= 4Jo t erfc ^)] ^ 

is found to be 

•-f(-f) 

Here is the work. Integration by parts is necessary. Define I = / 0 °° [erfc(?y )] 2 dr), such 
that e = jl. Let u = [erfc ^)] 2 and dv = dr). Then du — erfc(r;) e -7?2 ) dr) and 


v — r). Thus 


Since erfc(oo) = 0 , 


I — i) [erfc(77)]' 


OO 4 f°° __ 2 

\— 7= / V e v &k(r))dr) 

7T JO 


0 


4 r°° 2 

?= / 17 e 77 eiic(r))di) 

7T JO 




2 r°° 2 

= — 7 = / 2rje v eiicMdi) 

y/TT JO 

Integrate by parts again. Let w = erfc( 77 ) and dv = 2i) e~^di). Then du = — dr) 


and v = — e n . Hence, 


7 = ^[-e-’”erfc(,)|”--l^e-^ 


2 


-I fe-^dr) 

y/l r Jo 

Use the transformation w = \f2r). Then dw — \j2dr\. Substituting, 

2 1 r°° 2 

- — -= / e~ w dw 
\Ar v 2 Jo 


I = 


V* 

2 

V5F 


' \/2 V5F] 2 / _ V2\ 

V^F 2 J yfH [ L 2) 
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Therefore, 



f(-f) 
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APPENDIX C. VORTICITY FLUX 


The vorticity equation in the wall boundary layer is 


U U x + V U)y = oj 


XX 


This can be rewritten as (U ui) x + (vu)y = v xx using continuity. Now, integrate this 
with respect to x, obtaining 

d r°° 


o° d f°° 

U u> + / vujdx = u a 

x=o dY Jo 


oo 

x=0 


(C.l) 


As x —► oo, u —* 0 (and u x —*■ 0), and when x = 0, U = 0, so that 

d r°° 


d r 00 

—— / vudx = —oj x 

dY Jo *=o 


(C.2) 


The vorticity can also be written as u> = v x —Uy - Prom the continuity equation, 


U = - 



Vy dx 


which yields 

O ) = V X + — J Q Vy dx = V x + j Q Vyy dx. (C.3) 

A scahng analysis is performed on the vorticity terms. In the wall boundary 
layer, x ~ A ~ M~ l P~ x x, y ~ L ~ M~ l P~ 2 Y, u~P 2 U,v~Pv, and cD ~ MP 2 uj. 
Substituting into u — v x — Uy yields MP 2 u> = MP 2 v x + MP 4 Uy. The last term 
on the right is much smaller than the first term on the right (since Pel), giving 
u = v x . Thus, the vorticity u can be characterized by v x , ignoring the /q v yy dx term 
as small. This result is consistent with the numerical data of Canright [Ref. 1], as 
shown in Figure 17 (the vorticity plot is the solid line, the velocity derivative plot 
is the dotted line). As can be seen, the two plots are virtually identical for given 
horizontal slices of the vorticity and vertical velocity flux (in this case, at a depth of 
U = 0.0584). 
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Vorticity vs Velocity Flux 



Figure 17. Typical Vorticity and Vertical Velocity Flux Profiles at the same Depth 
from the Free Surface 

Back to the equation for the vorticity flux, 



= 0 (C.4) 

This shows that the vorticity flux at the wall is zero! Very unexpected. Sub¬ 
stituting in the expression for v and v x and then evaluating the integral does in¬ 
deed yield 0. Perhaps a higher order expansion for the velocity is needed, such as 
v « Pvi + P 2 v 2 -, to resolve the diffusion of vorticity from the wall. 
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APPENDIX D. TIMMAN’S METHOD ALONG 

THE WALL 


The velocity profile approximation for the wall boundary layer is a little more 
complicated than the surface profile approximation. The no-slip condition at the wall 
requires that the velocity is zero when x = 0. Also, the core velocity goes to 0, so 
the boundary condition as x —* oo requires a zero velocity. Between these two values 
of x, the velocity increases to the point in the boundary layer where the positive 
vorticity goes to zero, then the velocity decays through the boundary layer where the 
vorticity is negative (that vorticity which is convected down from the surface) until 
the velocity decays to zero. 

The following approximate velocity profile is assumed: 

v(x) = f(r}) = - f e~ t2 (a 4- ct 2 ) dt + e -7?2 (6 + drj 2 ), (D.l) 

Jr) 

where f(rj) = v, the vertical (dominant) component of the wall velocity, and 77 = 
where A is the thickness of the wall boundary layer. The boundary conditions are: 
at x = 0,77 = 0, and f(rj) = 0, due to the no-slip condition. As x —> 00, 77 —> 00, and 
/(» l) -* 0- 

The determination of the coefficients a, b, c, and d must be made. All four 
coefficients are functions of Y, the vertical direction. Start with the momentum 
direction in the Y-direction, 

Uv x + VVy = v xx (D.2) 

Integrate this equation with respect to x, from x = 0 at the wall to x = 00, outside 
the vertical boundary layer. This yields 


r°o I00 

/ (Uv x + VVy) dx = v x \ 

Jo 10 

Substitute the shear stress at the wall, r 0 , as r 0 = (j^J _ q = v x . Thus, 

roc 

— To = / ( Uv x + Wy)dx 

JO 


(D.3) 


(D.4) 
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(Uv) 


+ 



dx 


d rcc 

dY 


r oo 

/ »= 

Jo 


dx 


Define the momentum thickness to be A2, where 

f°° 

A 2 = / v 2 dx. 

J x=0 


Then 


roc roo 

A 2 = / v 2 dx = f 2 dx. 

Jx —0 Jx—O 


Since rj = ■%-, dx = A drj, and 


roo 

A 2 = A / v 2 dr}, 

Jx —0 


with A being the bmmdary-layer thickness. Therefore, 

dA 2 


-t 0 = 




(D. 5 ) 


(D.6) 


(D. 7 ) 


(D.8) 


The derivative conditions on the function f(rj) evaluated at the wall are as 
follows: 

'<"> - s -x( a+ i) 

m = « 

/"(O) = 2 (d-b) 
f"( 0) = 2(c — a) 

/ iv (0) = 12(6-2d) 


Since / = v, f( 0 ) = 0 , as the velocity at the wall is zero. Since there is no pressure 
gradient term in the momentum equation, the pressure is constant, and /"(0) = 0. 
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In addition, /"'( 0 ) = 0 , from the given boundary conditions [Ref. 31 ]. These last 
two conditions yield a = c and b = d. The first condition gives a = 0ne 

more condition is needed to solve for the four coefficients. The value of the shear 
stress at the wall is unknown, so /'( 0 ) is not specified, and the expression for / iv ( 0 ) is 
needed. Taking the fourth derivative of / with respect to rj yields f w (r}) = A 4 v x v x y. 
Since = Av x , then a = Av x , evaluated at x = 0 , and the last condition yields 
f"(r,)= A 3 o^(f). 

Equating the various derivative conditions yields two nonlinear equations in¬ 
volving the boundary-layer thickness A and the coefficient a. These are 

-T„ = (A,) 

and 

- i26 = A3a ^(£) 

Evaluating the integral associated with A 2 numerically, the first equation can be 
worked to the form 


-To = 



-Vx 



a 

A 


d 

dY 


(r 2A ) 


— (a 2 Ay + 2dA Qy'j 


£ 

A 


— aAy -f- 2Au y 


while the second equation can be worked to a similar form 


-126 = —= 




(D.9) 
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9 Vtt 

"A" 


aAy — A<2y 


(D. 10 ) 


These equations may be further reduced to 


ay = “ 3 v^) 

(D-ll) 

and 


Ay = ——(1 + 6x/tt) 

Z\a 

(D. 12 ) 


Initial conditions need to be specified on the displacement thickness (call it 
the mass flux) and on the flux of vorticity. The mass flux and the vorticity flux must 
be matched to the free surface boundary layer, so more than one free parameter may 
be needed. 

If d = 0 , then from the initial conditions, 6 = 0 , which further means that 
a = — f. The velocity profile is then 

roo 0 

/(?7) = v = a / e _t (2 1 2 — 1 ) dtt. (D. 13 ) 

Jr) 

T his approach proved to be more difficult than expected and was abandoned 
in favor of the modification to Glauert’s method. 
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APPENDIX E. NUMERICAL CODES 


. May YOU CONVERGE WITHOUT DELAY. - Dan Zinder 


The computer programs listed here axe supplied on an “as is” basis, with no 
warrantees of any kind. The author bears no responsibility for any consequences of 
using this program. 
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/* this file: fluid.h - header definitions for cold corner */ 

#include <stdio.h> 

#include <math.h> 

#define maxx 3001 
#define maixxvec (maxx - 2) 

#define Ma 10000 
#define Pr 0.01 

#define eps (sqrt(3.14159265)*(1.0-sqrt(2.0)/2.0)/2.0) 

#ifdef MAIN 

int i, j, k, nx, ny; 
int size; 

double hx, hy, dt; 

double x[maxx], y[maxx], b[maxx], xx[maxx], f [maxx] ; 
double xconl, xcon2, yconl, ycon2; 

double cl [maxx] [maxx], c2, c3[maxx] [maxx] , c4[maxx] [maxx] , c5, 
c6[maxx] [maxx] ; 

double dl [maxx] [maxx] , d2, d3 [maxx] [maxx] , d4 [maxx] [maxx] , d5, 
d6 [maxx] [maxx] ; 

double t [maxx] [maxx] , temp [maxx] [maxx] , tstar [maxx] [maxx] ; 

double hofx[maxx], integrand [maxx] , gmfun[maxx]; 

double gpsi[maxx] [maxx] , jpsi[maxx] [maxx] , psi[maxx] [maxx] ; 

double wpsi [maxx] [maxx] , spsi [maxx] [maxx] , unifpsi [maixx] [maxx] ; 

double u [maxx] [maxx] , v [maxx] [maixx] ; 

double surfu, massflux, F, eta, g; 

double delta [maixx] , seta [maxx] ; 

double tempa[maxxvec] , tempb [maxxvec] , tempo [maxxvec] ; 
double tstarn[maxx] , Tx[maxx], Txx[maxx] ; 
double A [3] [maxx] , C [3] [maxx] ; 

#undef MAIN 
#else 

extern int i, j, k, nx, ny; 

extern int size; 

extern double hx, hy, dt; 

extern double x [] , y [] , b [] , xx [] , f [] ; 

extern double xconl, xcon2, yconl, ycon2; 

extern double cl [] [maxx], c2, c3 [] [maxx] , c4[][maxx], c5, c6[] [maixx]; 
extern double dl [] [maixx] , d2, d3 [] [maxx] , d4 [] [maxx] , d5, d6 [] [maxx] ; 
extern double t[][maxx], temp [] [maxx] , tstar [] [maxx] ; 
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extern double hofx[] , integrand [] , grnfun[] ; 

extern double gpsi [] [maxx] , jpsi [] [maxx] , psi [] [maxx] ; 

extern double wpsi [] [maxx] , spsi [] [maxx], unifpsi [] [maxx]; 

extern double u [] [maxx] , v [] [maxx] ; 

extern double surfu, massflux, F, eta, g; 

extern double delta [] , seta[] ; 

extern double tempa[] , tempb[], tempo []; 

extern double tstam[] , Tx[] , Txx[]; 

extern double A [] [maxx] , C [] [maxx] ; 

#endif 
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/* this file: main.c - main program for cold corner */ 

#define MAIN 
#include "fluid.h" 

#define LINESIZE 128 

#define eps (sqrt(3.14159265)*(1.0-sqrt(2.0)/2.0)/2.0) 

main(argc,argv) 

int argc; char *argv[]; 

•C 

double adi(), tol=0.001, norm; 
int n, m, ninner=6, nouter=6; 

if ( argc > 1 ) sscanf ( argv[l] , '“/.d", feninner ); 

if ( argc > 2 ) sscanf ( argv [2] , "’/.d", ftnouter ); 

if ( argc > 3 ) sscanf ( argv[3], "/ilf", &tol ); 

n = 10; 
adiinit(); 

for (m = 1; n > 2 && m <= nouter; m++){ 
tx(); 

strength (); 
jet(); 
green(); 

for (i = 0; i <= nx; i++){ 
for (j = 0; j <= ny; j++){ 

psi [i] [j] = jpsi [i] [j] + gpsi [i] [j] ; 

> 

> 

for (i = 0; i <= nx; i++){ 
for (j = 0; j <= ny; j++){ 
if (i == 0 & j != 0){ 

u[i] [j] = Cpsi[i] [j+1] - psi[i] [j-1]) / ((2.0)*(y[j+l] -y[j])); 

v[i] [j] = - (psi[i+l] [j] - psi[i] [j]) / (x[i+l] - x[i]);> 

else if (j == 0 & i != 0){ 

u[i] [j] = (psi [i] [j+1] - psi[i] [j]) / (y[j+l] - yCj3); 
v[i][j] = - (psi[i+l][j] - psi[i-l][j]) / ((2.0)*(x[i+l] - x[i]));> 
else if (i == nx){ 

u[i] [j] = (psi[i] [j+1] - psi[i] [j-1]) / (2.0*(y[j+l] - y[j])); 

v[i] [j] = - (psi [i] [j] - psi [i-1] [j]) / (x[i] - x[i-l]);> 
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else if (j == ny){ 

u[i][j] = (psi [i] [j] - psi [i] [j-1]) / (y[j] - y[j-l]); 

v[i] [j] = - (psi [i+1] [j] - psi [i-1] [j]) / (2.0*(x[i+l] - x[i])); 

u [nx] [ny] = v [nx] [ny] = 0.0; > 

else{ 

u[i][j] = (psi[i] [j+1] - psi[i] [j-1]) / ((2.0)*(y[j+l] - y[j])); 
v[i] [j] = - (psi [i+1] [j] - psi [i-1] [j]) / ((2.0)*(x[i+l] - x[i]));> 
/* printf ( "7.1.2f 7.1.2f %1.5f 7.1.5f\n", x[i], y[j], u[i][j], v[i][j]>; */ 

> 

> 

norm = 1.0; 

for (n = 1; norm > tol && n < ninner; n++){ 
norm = adi(); 

/* printf ("Iteration number 7.d, norm: 7.f \n" ,n,norm) ; */ 

y 

fprintf (stderr, "Iteration = 7.d\n", m) ; 

> 

for (i = 0; i <= nx; i++){ 
for (j = 0; j <= ny; j++){ 

printf ("74.3f %1.3f 7.1.6f\n\ x[i] , y[j], temp[i][j]); 

> 

> 
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/* this file: adiinit.c - initialization for cold corner */ 
#include "fluid.h" 
adiinitO { 


nx = ny = maxx - 1; 
hx = hy = 3.0 / nx; 
dt = 0.4 * hx; 

xconl = dt / (hx * hx); 
xcon2 = dt / (4.0 * hx); 
yconl = dt / (hy * hy); 
ycon2 = dt / (4.0* hy); 

/* printf("7,1.4f 7,1.4f 7.1.4 7,1.4f\n", xconl, xcon2, yconl, ycon2); */ 

c2 = 1.0 + xconl; 

c5 = 1.0 - yconl; 

d2 = 1.0 + yconl; 

d5 = 1.0 - xconl; 

for (i = 0; i <= nx ; i++){ 
x[i] = y[i] = i * hx; 
f[i] = 0.0; 

> 

/* printf(" initial\n"); 

printf ("x[i] y[j] temp[i] [j]\n") ; */ 
for (i = 0; i <= nx ; i++H 
for (j = 0; j <= ny ; j++){ 
if (i == 0) 

temp[i] [j] = tstar[i] [j] = -1.0; 
else if (i == nx) 

temp[i] [j] = tstar[i][j] = 0.0; 
else if (j == 0) 

temp[i] [j] = -exp(-10.0 * x[i]); 
else 

temp[i][j] = tstar[i][j] = 0.0; 

/* printf ("7,1.3f 7,1.3f 7,1.4f\n", x[i], y[j] , temp[i][j]); */ 

> 

> 

} 
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/* this file:adiinit2.c - initialization for adi solver at fine grid. */ 


#include "fluid.h" 
#define LINESIZE 128 

adiinit() { 


nx = ny = maxx - 1; 
hx = hy = 3.0 / nx; 
dt = 0.4 * hx; 

xconl = dt / (hx * hx); 
xcon2 = dt / (4.0* hx); 
yconl = dt / (hy * hy); 
ycon2 = dt / (4.0 * hy); 

/* printf ("‘/,1.4f 7,1.4f 7.1.4 7.1.4f\n", xconl, xcon2, yconl, ycon2); */ 

c2 = 1.0 + xconl; 
c5 = 1.0 - yconl; 
d2 = 1.0 + yconl; 
d5 = 1.0 - xconl; 

for (i = 0; i <= nx ; i++){ 
x[i] = y[i] = i * hx; 
f[i] = 0.0; 

} 

for (i = 0; i <= nx; i++){ 
for(j = 0; j <= ny; j++){ 

scanf ("%*lf7.*lf7.1f" , &(temp[i] [j])) ; /* Read in converged data */ 

tstar[i][j] =temp[i][j]; 

> 

> 

for (j = 0; j <= ny ; j++){ 

temp[0] [j] = tstar[0] [j] = -1.0; 
temp[nx] [j] = tstar[nx] [j] = 0.0; 

> 

> 
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/* this file: tx.c This file calculates the thermal gradients, 

Tx and Txx, and the heat flux. */ 

#include "fluid.h" 

tx() 

for (i = 0; i <= nx; i++M 
if (i == 0) 

Tx[i] = (temp [i+1] [0]-temp [i] [0]) / (x[i+l]-x[i]); 
else if (i == nx) 

Tx[i] = (temp [i] CO]-temp [i-1] [0]) / (x[i]-x[i-l]) ; 
else 

Tx[i] = (temp[i+1] [0]-temp[i-1] [0] ) / (2.0* (x[i+l]-x[i]>) ; 
/* printf ("7,1.2f %1.4f 7.1.4f\n", x[i] , temp[i][0], Tx[i]); */ 

> 

for (i = 0; i <= nx; i++){ 
if (i == 0) 

Txx [i] = (Tx[i+1] - Tx[i]) / (x[i+l]-x[i]); 
else if (i == nx) 

Txx[i] =0; 
else 

Txx[i] = (temp [i+1] [0]-2*temp[i] [0]+temp [i-1] [0] ) / 
pow((x[i+l]-x[i]),2.0); 

/* printf ("%1.2f %1.4f /Cl.4f\n", x[i] , temp[i][0], Txx[i]); */ 

> 

for (i = 0; i <= nx; i++){ 

f[i] = (-1.0/2.0)* pow((fabs(-temp[i] [0]*Tx[i])/eps),(2.0/3.0)); 

> 
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/* this file: strength2.c This is a file to determine the 
characteristic velocity from the surface boundary layer into the corner. */ 

#include "fluid.h" 

strength() 

{ 

F = fabs((16.0/5.0)*(sqrt(3.14159265))*(pow(Tx[0]/(eps*2.0) , 1.0/3.0))) ; 

/* printf ("F = 7.1.6f\n", F); */ 

> 
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/* this file: jet.c This file determines the velocity from the 
potential flow due to a wall jet in the quarter-plane. */ 

#include "fluid.h" 

jetO 

double theta, con; 
con = (5.0/2.0) * F; 

/* printf("x[i] yCj] jpsi [i] [j]\n") ; */ 
for (i = 0; i <= nx; i++){ 
for(j = 0; j <= ny; j++){ 
theta = atan2(x[i] ,y [j] ) ; 
if (j == 0) 

jpsi[i] [j] =0.0; 
else if (i == 0) 

jpsi[i][j] = -2.0 * pow(con, 0.25) * powCy[j],0.25); 
else 

jpsi[i] [j] = -2.0 * pow(con, 0.25) * pow((x[i]*x[i]+y[j]*y[j]), 
0.125) * (cosCtheta/4.0) - (1.0+sqrt(2.0))*sin(theta/4.0)); 

/* printf ('7,1.2f 7.1.2f 7,1.4f\n", x[i], y[j] , jpsi[i] [j]) ; */ 

> 

> 

} 
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/* this file: green.c This file solves Laplace’s equation in core. */ 


#include "fluid.h" 

greenO 

{ 

double piece; 

for (i = 0; i <= nx; i++){ 
if (i < nx) 

hofx [i] = ( pow( ((-temp[i] [0])*(-temp[i][0])/(eps*eps*Tx [i])) , 

1.0/3.0) ) / (-2.0); 
else 

hofx[i] =0.0; 

/* printf("%1.2f */,1.12f\n", x[i], hofx[i]); */ 

} 

for (i = 1; i <= nx; i++){ 
f°r(j = 1; j <= ny; j++){ 
piece = 0.0; 
integrand[0] = 0.0; 
for (k = 1; k <= nx; k++){ 

gmfun[k] = ( y[j]/((x[k]+x[i])*(x[k]+x[i]) + y[j]*y[j]) - 
yCj]/((x[k]-x[i])*(x[k]-x[i]) + y[j]*y[j]) )/( 3 . 14159265) ; 
integrand [k] = (- grnfunCk] * hofx[k]); 

piece += C(x[k]-x[k-l] ) / 2.0) * (integrand[k]+integrand[k-l] ) ; 

/* printf ("*/,1.4f\n", integrand [k] ) ; */ 

> 

gpsi [i][j] = piece + (hofx[0] * (2.0/3.14159265) * atan2(y[j],x[i])); 

> 

> 

for (i = 0; i <= nx; i++){ 
gpsi[0] [i] = hofx[0] ; 
gpsi[i] [0] = hofx[i] ; 

> 

> 
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/* this file: adi.c - solve the heat equation for the potential flow 
due to a wall jet in the cold corner (quarter plane). 

2 2 


d T 

d T 

d T 

d T 

d T 

M(- + 

u- + 

v-) 

=- + 

— 

d t 

d x 

d y 

2 

2 


d x d y 

with BC: 

at y = 0: Ty = 0, uy = Tx, v = 0 

at x = 0: T = -1, u = v = 0 

using the symmetric ADI method, which uses two subprograms for computing 
the L U decomposition of a tridiagonal matrix and then for solving 
L U x = b. A system of equations then result of the form (a two-step 
process): 


* n 

(Step 1) A * T = B * T 

x,y x,y 

n+1 * 

(Step 2) C * T = D * T 

x,y x,y 


where A, B, C, and D are all tridiagonal matrices. 

*/ 

#include "fluid.h" 
double adiO 

int i, j, m, n; 
double max, diff; 

/* Next, bring in the values of u (psi_y) and v (- psi_x) from the 
potential flow due to a wall jet. Input the velocities from file 
main.c. */ 

/* STEP 1: Solve the tridiagonal system 
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* * * n n n 

cl T + c2 T + c3 T = c4 T + c5 T + c6 T 

x+l,y x,y x-l,y x,y+l x,y x,y-l 

where the coefficients cl through c6 are determined by */ 

for (i = 0; i <= nx ; i++){ 
for (j = 0; j <= ny ; j++){ 

cl[i][j] = (xcon2 * u[i][j]) - (xconl / 2.0); 
c3[i][j] = (-xcon2 * u[i][j]) - (xconl / 2.0); 
c4[i] [j] = (-ycon2 * v[i][j]) + (yconl / 2.0); 
c6[i][j] = (ycon2 * v[i][j]) + (yconl / 2.0); 

/* printf("*/,d 7.d 7.1 Af %1.2f %1.4f 7.1.4f 7.1.2f 7.1.4f\n\ i, j, 
cl [i] [j] , c2, c3 [i] [j] , c4[i][j], c5, c6[i][j]); */ 

> 

> 

/* c2 and c5 are defined in initialization file. */ 

/* Start with the corner. */ 

j = 0; 
i - 1; 

tempa[i-l] = yconl * temp[i] [j+1] + c5 * temp[i] [j] - (yconl * 
hy * f[i]) - c3[i][j] * temp[i-l] [j]; 

A[l] [i-1] = c2; 

A[2] [i-1] = cl [i] [j] ; 

/* printf("7.d 7.d %1.4f\n", i, j, tempa[i-l]); */ 

/* Next, go along the surface to the other corner. */ 
for (i=2;i<nx-l; i++){ 

tempa[i-l] = yconl * temp[i] [j+1] + c5 * temp[i] [j] - (yconl * 
hy * f [i]); 

A[0] [i-1] = c3[i] [j] ; 

A[l] [i-1] = c2; 

A [2] [i-1] = cl[i] [j] ; 

/* printf("7.d 7.d 7,1.4f\n", i, j, tempa[i-l]); */ 

> 

i = nx - 1; 

tempa[i-l] = yconl * temp[i][j+l] + c5 * temp[i] [j] - (yconl * 
hy * f[i]) - cl [i] [j] * temp[i+l] [j] ; 
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A[0] [i-1] = c3[i] [j] ; 

A[1][i-1] = c2; 

/* printf("7od 7.d %1.4f\n", i, j, tempa[i-l]); */ 

/* for (m = 0; m < 3; m++H 

for (n = 0; n < ny - 1; n++){ 

printf("7,d 7.d 7,1.4f\n", m, n, A[m][n]); 

> 

} */ 

size = maxxvec; 

/* Use tridiag.c to solve A T* = Tn, where A = is tridiagonal, 

Tn is the known temperature (tempa), and T* is the step temperature. */ 

tridiag(A, tempa, size); 

for (i = 1; i <= size ; i++){ 
tstar[i][j] = xx[i-1] ; 

/* printf ("7.d 7.d 7.1.4f\n", i, j, tstar[i] [j]); */ 

} 

/* This is the solution to the surface row. The tstar matrix is 
updated with the new surface temperature. Now, work on the inside 
of the regime. Continue to update the tstar matrix as each row is 
solved. */ 

for (j = 1; j < ny; j++M 
i = l; 

tempb[i-1] = c4[i] [j] * temp[i] [j+1] + c5 * temp[i][j] + c6[i] [j] 

* temp[i] [j-1] - c3[i] [j] * temp[i-l] [j] ; 

A[l] [i-1] = c2; 

A [2] [i-1] = cl [i] [j] ; 

/* printf ("7.d ‘/,d 7.1.4f\n", i, j, tempb [i-1]); */ 
for (i=2;i<nx-l; i++M 

tempb [i-1] = c4[i] [j] * temp [i] [j+1] + c5 * temp[i] [j] + c6[i] [j] 

* temp[i] [j-1] ; 

A[0] [i-1] = c3[i] [j] ; 

A[l] [i-1] = c2; 

A [2] [i-1] = cl [i] [j] ; 

> 
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i = nx - 1; 

tempb[i-l] = c4[i] [j] * temp[i] [j+1] + c5 * temp[i][j] + c6[i] [j] 

* temp[i][j-l] - cl[i][j] * temp[i+l] [j] ; 

A [0] [i-1] = c3[i] [j] ; 

A Cl][i-1] = c2; 

/* printf("7,d %d 7.1.4f\n", i, j, tempb[i-l]); */ 

/* Use crout(A) and croutslvCL, U, tempb) to solve for row j. Then 
update the temp matrix. */ 

tridiag(A, tempb, size); 

for (i = 1; i <= size ; i++){ 
tstar[i][j] = xx[i-l]; 

/* printf ("%d 7.d 7,1.4f\n", i, j, tstar [i] [j] ) ; */ 

> 

> 

/* Finally, solve along the bottom row. */ 

j = ny; 
i = 1; 

tempo [i-1] = c5 * temp[i] [j] + yconl * temp [i] [j-1] - c3[i] [j] 

* temp [i-1] [j] ; 

A[l] [i-1] = c2; 

A [2] [i-1] = cl[i] [j] ; 

/* printf ("*/,d */,d */,1.4f\n", i, j, tempo [i-1]); */ 

for (i = 2; i < nx - 1; i++M 

tempo [i-1] = c5 * temp[i] [j] + yconl * temp[i] [j-1] ; 

A[0] [i-1] = c3[i] [j] ; 

A[l] [i-1] = c2; 

A [2] [i-1] = cl [i] [j] ; 

/* printf ("7,d °/,d 7,1.4f\n", i, j, tempo [i-1]); */ 

y 

i = nx - 1; 

tempc[i-l] = c5 * temp[i] [j] + yconl * temp[i] [j-1] - cl[i][j] 

* temp [i+1] [j] ; 

A[0] [i-1] = c3[i] [j] ; 

A[l] [i-1] = c2; 

/* printf ("7,d 7.d 7.1.4f\n", i, j, tempo [i-1]); */ 
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for (m = 0; m < 3; m++){ 

for (n = 0; n < ny - 1; n++){ 

/* printf ("7.d %d %1.4f\n", m, n, A[m] [n] ) ; */ 

> 

> 

for (i = 0; i < nx - 1; i++){ 

/* printf ('7.1.4f\n\ tempo [i]); */ 

> 

tridiag(A, tempo, size); 

for (i = 1; i <= size ; i++){ 
tstar[i] [j] = xx[i-l] ; 

/* printf ("7.d 7.d %1.4f\n", i, j, tstar [i] [j]) ; */ 

y 

/* Print the results of the first step of tstar. */ 

/* printf C"x[i] y[j] tstar[i] [j]\n") ; */ 
for (i = 0; i <= nx; i++){ 
for (j - 0; j <= ny; j++){ 

/* printf ("*/,1.2f 7.1.2f 7.1.4f\n", x[i], yCj] , tstar [i] [j]) ; */ 

} 

} 

/* STEP 2: Repeat in the other direction. Now solve 

n+1 n+1 n+1 * * * 

dl T + d2 T + d3 T = d4 T + d5 T + d6 T 

x,y+l x,y x,y-l x+l,y x,y x-l,y 

where the coefficients dl through d6 are determined by */ 

for (i = 0; i <= nx ; i++){ 
for (j = 0; j<= ny ; j++){ 

dl[i] [j] = (ycon2 * v[i][j]) - (yconl / 2.0); 
d3[i][j] = (-ycon2 * v[i][j]) - (yconl / 2.0); 
d4[i][j] = (-xcon2 * u[i] [j]) + (xconl / 2.0); 
d6[i][j] = (xcon2 * u[i][j]) + (xconl / 2.0); 

/* printf("%d 7.d 7.1.4f 7.1-4f 7.1.4f 7.1-4f 7.1.4f 7.1.4f\n", i, j, 
dl [i] [j] , d2, d3[i][j], d4[i][j], d5, d6[i][j]); */ 

} 
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> 

/* d2 and d5 are defined in initialization file. */ 

/* Only one loop is needed, as the equations axe the same for all 
i columns. */ 

/* printf C"x[i] y[j] temp[i] [j]\n") ; */ 
size = maxx; 
max = 0.; 

for (i = 1; i < nx; i++){ 

j = 0; 

tstarn[j] = d4[i][j] * tstar[i+l] [j] + d5 * tstar[i] [j] + d6[i][j] 

* tstar[i-l] [j] - (yconl * hy * f[i]); 

C[l] [0] = d2; 

C [2] [0] = -yconl; 

for (j = 1; j < ny; j++){ 

tstarn[j] = d4[i] [j] * tstar[i+l] [j] + d5 * tstar[i] [j] + d6[i] [j] 

* tstar [i-1] [j] ; 

CEO] Cj] = d3[i] [j] ; 

C[l] [j] = d2; 

C[2] [j] = dl [i] [j] ; 

> 

j = ny; 

tstarn[j] = d4[i] [j] * tstar[i+1] [j] + d5 * tstar[i] [j] + d6[i][j] 

* tstar [i-1] [j] ; 

C[0][ny] = -yconl; 

C[l] [ny] = d2; 

tridiag(C, tstarn, size); 

for (j - 0; j < size ; j++){ 

diff = fabsC temp[i][j] - xx[j]); 
if (diff > max) { 
max = diff; 

/* printf ("new max [*/,d,7,d] = 7,f\n" ,i, j ,max) ; */ 

} 

temp [i] [j] = xx[j] ; 

> 

> 

/* Print the results for the second step and begin the next iteration 
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with the new temp. */ 


for (i = 0; i <= nx; i++){ 
for (j = 0; j <= ny; j++){ 

/* printf("*/,1.2f */,1.2f */,1.4f\n", x[i], y[j], temp[i][j]>; */ 

> 

> 


return(max); 




/* this file: tridiag.c Solve tridiagonal system A x = b. */ 

#include "fluid.h" 

tridiag(t, b, size) 

double t[3][maxx], b[maxx]; 
int size; 

{ 

double r; 

for (j = 1; j< size; j++){ 
r = t[0][j] / t [1] [j—1] ; 
t Cl] [j] -= r * t [2] [j-1] ; 
b[j] -= r * b[j-l] ; 

> 

b[size-1] /= t[l] [size-1] ; 

/* printf("7,d %1.4f\n", size-1, b[size-l]); */ 

for (j = size-2; j >= 0; j—){ 

b[j] = (b[j] - (t[2] [j] * b[j+l])) / t[1] [j] ; 

/* printf("7,d 7,1.4f\n", j, b[j]); */ 

> 

/* Solution is stored in xx[]. */ 

for (j = 0; j< size; j++){ 
xx[j] = b[j] ; 

> 

> 
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/* this file: surface.c This is a file to determine the stream 

function from the surface boundary layer. */ 


#include "fluid.h" 
surface() 

for (i = 0; i < nx; i++){ 

delta[i] = pow(eps,-l.0/3.0) * pow(-temp[i][0] ,1.0/3.0) * 
pow(Tx[i],1.0/3.0); 

seta[i] = - pow(eps,-2.0/3.0) * pow(-temp[i][0],2.0/3.0) * 
pow(Tx[i],-l.0/3.0); 

for (j = 0; j <= ny; j++){ 

spsi[i][j] = seta[i] * ( (1.0/2.0) * (1.0 - 
exp(- pow(y[j]/(Pr*delta[i]),2.0))) + ((sqrt(3.14159265)/2.0) 

* (y [j]/(Pr*delta[i])) * erfc(y[j]/(Pr*delta[i]))) ); 

/* printf("%1.3f %1.3f %1.6f\n", x[i], y[j], spsi[i][j]); */ 

y 

y 

y 
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/* this file: wall.c This is a file to determine the stream 
function in the wall boundary layer. */ 

#include "fluid.h" 

wall() 

double coni, con2, psicon, Delta, feta; 
double getgO; 

coni = (5.0/2.0) * F; 

con2 = 2.0 / pow(coni,1.0/4.0); 

for (j = 1; j <= ny; j++){ 

Delta = con2 * pow(y[j],3.0/4.0); 
psicon = - 2.0 * pow((conl*y[j]),1.0/4.0); 

for (i = 0; i <= nx; i++){ 
eta = x[i] / (Delta * Pr); 

/* printf("eta = */.1.3f\n", eta); */ 
g = getg(eta); 
feta = g * g; 

/* printf ("feta = */.1.6f\n", feta); */ 
wpsi[i][j] = psicon * feta; 

/* printf ("7.1.3f %1.3f 7.1.6f\n", x[i] , y[j] , wpsi[i][j]); */ 

> 

> 

> 
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/* This file: getg.c - This is a file to determine eta using 
Newton's Method, to be used in wall.c */ 

#include "fluid.h" 

double getg(eta) 

double eta; 

double gO, gprime, etaO, tol; 

tol = 0.000001; 

if (eta > 7.8) return(l.O); 

if (eta < exp(1.0)) 
gO = eta / 3.0; 
else 

{gO = 1.0 - sqrt(3.0 / exp(2.0 * eta - 3.14159265 / sqrt(3.0)));} 
/* printfC'gO = */,1.6f\n", gO); */ 


g = gO; 
do{ 

gO = g; 

gprime = (1 - g * g * g)/ 3.0; 

etaO = Iog((sqrt(1.0+g0+g0*g0))/(1.0-g0)) + (sqrt(3.0)) * 
atan(((sqrt(3.0))/3.0)*(l.0+2.0*g0)) - (sqrt(3.0))*(3.14159265/6.0); 

g = gO + (eta - etaO) * gprime; 

/* printf ("g = */,1.6f\n", g); */ 

} while (fabs(g - gO) > tol); 
return(g); 


> 
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/* this file: main.c - main program for uniform stream function. */ 


#define MAIN 
#include "fluid.h" 

#define LINESIZE 128 

#define eps (sqrt(3.14159265)*(1.0-sqrt(2.0)/2.0)/2.0) 
main() 

-c 

psiinit(); 

txO; 

strength(); 

jet(); 
greenC); 
wallO; 
surface(); 

for (i = 0; i <= nx; i++){ 
for (j = 0; j <= ny; j++){ 

unifpsi [i] [j] = jpsi[i][j] + gpsi[i][j] + wpsi[i][j] 

+ spsi [i] [j] - jpsi [0] [j] - gpsi [i] [0] ; 

printf ("y,1.3f */.1.3f */,1.6f %1.6f */.1.6f */,1.6f */,1.6f\n", 
x[i], y[j], jpsi [i] [j], gpsi[i][j], spsi[i][j], wpsi[i][j], 
unifpsi [i] [j]); 

> 

> 

> 
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/* this file: psiinit.c - initialization for uniform stream function. */ 


#include "fluid.h" 
#define LINESIZE 128 

psiinitO { 

FILE *file; 

char line[LINESIZE]; 

nx = ny = maxx - 1; 
hx = hy = 3.0 / nx; 


for (i = 0; i <= nx ; i++){ 
x[i] = y[i] = i * hx; 

> 

file = fopen("input.txt", "r"); 

for (i = 0; i <= nx; i++){ 
forCj = 0; j <= ny; j++){ 

fscanf (f ile, "°/.lf*/.lf7,lf\n" , &Cx[i]), &(y[j]), &(temp[i] [j] )) ; 

> 

> 

/* printf C'"/,1.3f 7,1.3f %1.6f\n", x[l], y[2], temp [1] [2]); */ 

} 
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/* this file:unitemp.c - Calculates the uniform temperature solution 
by matching Tinner + Touter - Tmatch. */ 

#define MAIN 
#include "fluid.h" 

#define LINESIZE 128 

#define eps (sqrt(3.14159265)*(1.0-sqrt(2.0)/2.0)/2.0) 

mainC) 

{ 

double unitemp [maxx] [maxx] , teta; 

nx = ny = maxx - 1; 
hx = 5.0 / nx; 

for (i = 0; i <= nx ; i++){ 
x[i] = y[i] = i * hx; 

> 

for (i = 0; i <= nx ; i++){ 
for(j = 0; j <= ny; j++){ 

scanf ("7.*lf 7.*lf 7.1f " , &(temp [i] [j] ) ) ; 

> 

> 

for (i = 0; i <= nx; i++){ 
if (i == 0) 

Tx[i] = (4.0*temp[l] [0]-3.0*temp[0] [0]-temp[2] [0] )/(2.0*x[i+1]) ; 
else if (i == nx) 

Tx[i] = (temp [nx] [0]-temp [nx-1] [0]) / (x[i]-x[i-l] ) ; 
else 

Tx[i] = (temp[i+1] [0]-temp[i—1] [0]) / (2.0* (x[i+l]-x[i] )) ; 

/* printf("7.1.2f 7.1.4f 7.1.4f\n\ x[i], temp[i][0], Tx[i]); */ 

> 

for (i = 0; i <= nx ; i++){ 
f°r(j = 0; j <= ny; j++){ 

teta = y[j] / (Pr * pow((-temp[i] [0])/(eps * Tx[i] * Tx[i]), 

1.0/3.0) ); 

/* printf("7.1.2f 7.1.2f 7,1.15f\n", x[i], y[j], teta); */ 

unitemp[i] [j] = temp[i] [j] + (Pr / (4.0*eps)) * temp[i][0] * 

( ((sqrt(3.14159265)) * (erfc(teta))) * (0.5 + teta * teta) - 
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(teta * exp(- teta * teta)) ); 

printf ("%1.2f %1.2f %1.5f\n", x[i], y[j], unitemp [i] [j]); 

> 

> 

} 
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